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13.  ABSTRACT  {Maximum  200  words) 


This  report  provides  a  description  of  the  Navy  Coastal  Ocean  Model  (NCOM)  Version  1.0.  The  model  has  a  free  surface  and  is 
based  on  the  primitive  equations  and  the  hydrostatic,  Boussinesq,  and  incompressible  approximations.  The  model  uses  an  Arakawa 
C  grid,  is  leapfrog  in  time  with  an  Asselin  filter  to  suppress  timesplitting,  and  uses  second-order  centered  spatial  finite  differences. 
The  propagation  of  surface  waves  and  vertical  diffusion  is  treated  implicitly.  A  choice  of  the  Mellor-Yamada  Level  2  or  Level  2.5 
turbulence  models  is  provided  for  the  parameterization  of  vertical  mixing. 

The  horizontal  grid  is  curvilinear.  The  vertical  grid  uses  sigma  coordinates  for  the  upper  layers  and  z-level  (constant  depth) 
coordinates  for  the  lower  layers,  and  the  depth  at  which  the  model  changes  from  sigma  to  z-level  coordinates  can  be  specified  by  the 
user.  The  combined  vertical  coordinate  system  provides  some  flexibility  in  setting  up  the  vertical  grid  and  easily  allows  comparisons 
to  be  made  between  simulations  conducted  with  sigma  and  z-level  coordinates.  The  inclusion  of  a  source  term  in  the  model  equations 
simplifies  input  of  river  and  runoff  inflows.  Some  limitations  of  the  model  are  discussed. 
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DESCRIPTION  OF  THE  NAVY  COASTAL  OCEAN  MODEL  VERSION  1.0 


1.  INTRODUCTION 

This  report  presents  a  description  of  the  Navy  Coastal  Ocean  Model  (NCOM)  Version  1.0.  The 
report  discusses  the  physics  and  numerics  of  the  model  and  some  of  its  limitations. 

The  goals  in  developing  NCOM  were  to  provide  an  ocean  model  that  would  include  the  best 
features  of  existing  ocean  models,  would  meet  the  Navy’s  needs  for  coastal  ocean  simulation  and 
prediction,  and  would  be  fully  compatible  with  and  could  be  fully  coupled  with  the  Navy’s  Coupled 
Ocean  and  Atmospheric  Mesoscale  Prediction  System  (COAMPS)  (Hodur  1997). 

The  first  of  these  goals  is  somewhat  elusive  in  that  new  ocean  model  features  are  continually 
being  developed  within  the  ocean  modeling  community.  The  approach  taken  for  the  initial  devel¬ 
opment  of  NCOM  1.0  is  to  make  use  of  fairly  well-established  ocean  modeling  techniques  that  have 
been  demonstrated  to  work  well  and  to  incorporate  improvements  and  additional  capabilities  into 
NCOM  when  they  are  determined  to  be  useful  or  needed.  Also,  it  may  not  be  possible  to  meet  all 
the  Navy’s  coastal  modeling  needs  with  a  single  model.  However,  the  approach  is  to  make  NCOM 
as  flexible  as  possible  without  incurring  a  significant  penalty  in  terms  of  efficiency. 

COAMPS  was  developed  by  the  Marine  Meteorology  Division  of  the  Naval  Research  Labora¬ 
tory  (NRL)  at  Monterey,  California,  and  has  a  very  specific  code  structure  that  provides  for  calls 
to  (i.e.,  coupling  between)  both  the  atmospheric  and  ocean  models  within  the  same  Fortran  pro¬ 
gram.  COAMPS  also  provides  for  an  arbitrary  number  of  levels  of  nesting  within  the  same  Fortran 
program.  This  nesting  capability  is  made  possible  by  using  dynamic  memory  allocation  with  array 
dimensions  specified  at  run  time  and  by  passing  model  variables  to  subroutines  through  subroutine 
argument  lists  rather  than  through  common  blocks.  This  allows  the  same  model  routines  to  calcu¬ 
late  the  different  nests.  Since  most  ocean  models  are  not  structured  in  this  way,  a  certain  amount 
of  recoding  is  required  to  fully  adapt  an  existing  ocean  model  to  the  COAMPS  code  structure. 

NCOM  Version  1.0  is  based  primarily  on  two  existing  ocean  models,  the  Princeton  Ocean 
Model  (POM)  and  the  Sigma/Z-level  Model  (SZM).  POM  was  developed  by  Alan  Blumberg  and 
George  Mellor  (Blumberg  and  Mellor  1983;  Blumberg  and  Mellor  1987).  POM  is  well-known  to 
anyone  familiar  with  ocean  models  and  has  attracted  a  wide  base  of  users  within  the  academic, 
civilian,  and  government  communities.  Table  1  lists  some  of  the  main  features  of  POM.  POM  is  a 
three-dimensional  (3-D),  primitive  equation,  baroclinic,  hydrostatic,  free  surface  model  and  uses  an 
orthogonal-curvilinear  horizontal  grid,  a  sigma  (i.e.,  bottom-following)  vertical  grid,  a  split-explicit 
treatment  of  the  free  surface,  Smagorinsky  horizontal  mixing,  and  the  Mellor-Yamada  Level  2.5 
(MYL2.5)  turbulence  model  for  vertical  mixing. 

Manuscript  approved  August  3,  1998. 
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Table  1  —  Comparison  of  Features  of  POM  and  SZM 


POM 


SZM 


Primitive  Equation 
Incompressible 
Hydrostatic 
Boussinesq 
Free  Surface 

Smagorinsky  Horizontal  Mixing 
Mellor-Yamada  Level  2.5  Vertical  Mixing 
Quadratic  Bottom  Drag 
Curvilinear  Horizontal  Grid 
Sigma  Coordinate  Vertical  Grid 
Leapfrog  in  Time  with  Asselin  Filter 
Second-Order,  Centered  Finite  Differences 
Flux  Conservative  Formulation 
Split-Explicit  Treatment  of  Free  Surface 
Implicit  Vertical  Mixing 


Primitive  Equation 
Incompressible 
Hydrostatic 
Boussinesq 
Free  Surface 

Grid- Cell  Re  Number  Horizontal  Mixing 
Mellor-Yamada  Level  2  Vertical  Mixing 
Quadratic  Bottom  Drag 
Cartesian  Horizontal  Grid 
Combined  Sigma/Z-level  Vertical  Grid 
Leapfrog  in  Time  with  Asselin  Filter 
Second-Order,  Centered  Finite  Differences 
Flux  Conservative  Formulation 
Implicit  Treatment  of  Free  Surface 
Implicit  Vertical  Mixing 


SZM  was  developed  in-house  at  NRL,  Stennis  Space  Center,  MS  (Martin  et  al.  1998).  SZM 
is  similar  in  many  ways  to  POM  (see  Table  1)  but  differs  from  POM  in  that  it  uses  a  Cartesian 
horizontal  grid  (the  grid  spacing  in  the  horizontal  is  constant),  a  combined  sigma/ 2-level  vertical 
grid  with  sigma  layers  near  the  surface  and  2-levels  (i.e,  fixed-depth  levels)  below  a  user-specified 
depth,  an  implicit  treatment  of  the  free  surface,  horizontal  eddy  coefficients  calculated  based  on 
a  maximum  grid-cell  Reynolds  number  criteria,  and  vertical  eddy  coefficients  calculated  using  the 
Mellor-Yamada  Level  2  (MYL2)  turbulence  closure  scheme. 

In  a  coastal  model  comparison  study  that  was  conducted  at  NRL  (Martin  et  al.  1998),  POM 
and  SZM  were  found  to  simulate  a  number  of  basic  physical  processes  that  are  important  in  the 
coastal  ocean  (advection,  mixing,  and  the  propagation  of  surface,  internal,  and  coastal  trapped 
waves)  fairly  well. 

NCOM  uses  those  features  of  POM  and  SZM  that  tend  to  provide  the  most  flexibility  and 
provides  a  choice  in  many  cases  where  selection  of  one  scheme  or  parameterization  over  another  is 
difficult  because  of  trade-offs  that  can  be  situation  dependent.  The  horizontal  grid  is  orthogonal- 
curvilinear  as  in  POM,  which  allows  adaptation  to  different  map  projections  and  also  allows  for  use 
of  a  limited  amount  of  grid  curvature  to  follow  a  slowly  curving  coastline  or  bathymetric  feature 
and/or  to  provide  increased  grid  resolution  in  certain  areas  of  the  domain.  The  sigma/ 2-level 
vertical  grid  from  SZM  is  used  to  offer  a  choice  of  sigma  layers  or  2-levels,  or  some  combination 
with  sigma  layers  in  the  shallow  water  and  2-levels  in  the  deeper  water.  A  choice  of  grid-cell 
Reynolds  number  or  Smagorinsky  horizontal  diffusion  and  MYL2  or  MYL2.5  vertical  diffusion  is 
provided. 

The  free  surface  is  treated  implicitly  as  in  SZM.  This  is  significantly  simpler  than  the  split- 
explicit  scheme  used  in  POM  and  is  more  consistent  in  the  sense  that  the  depth-averaged  equations 
are  the  almost  exact  vertical  integral  of  the  3-D  equations.  However,  the  implicit  scheme  is  not 
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as  accurate  for  the  propagation  of  surface  gravity  waves  because  of  the  much  larger  timestep  used 
to  propagate  these  waves  (for  some  comparisons,  see  Martin  et  al.  1998).  However,  an  implicit 
treatment  of  the  free  surface  has  been  used  in  a  number  of  coastal  models  (Leendertse  1989; 
Blumberg  1992;  Casulli  and  Cattani  1994;  Casulli  and  Cheng  1994).  An  implicit  treatment  of 
the  free  surface  has  usually  been  found  to  be  sufficiently  accurate  to  simulate  tides  and  wind  setup 
since  the  length  and  time  scales  of  these  processes  are  relatively  long.  It  was  originally  planned  to 
ofFer  a  choice  of  an  implicit  or  split-explicit  treatment  of  the  free  surface.  However,  there  are  some 
inconsistencies  between  the  use  of  a  split-explicit  scheme  and  a  z-level  type  of  grid  that  have  not 
yet  been  resolved.  The  option  of  a  split-explicit  treatment  of  the  free  surface  may  be  provided  at 
a  later  time. 

Some  additional  features  of  NCOM  are  (1)  a  source  term  in  all  the  equations  to  simplify  input 
of  rivers  and  coastal  runoff,  (2)  the  option  of  including  forcing  by  the  surface  air  pressure  and  the 
tidal  potential,  (3)  a  multicomponent  scalar  variable  array  that  allows  additional  scalar  fields  to 
be  easily  added  to  the  model  (e.g.,  for  optical  or  biological  submodels),  (4)  shrink-wrapping  to 
eliminate  calculations  over  land  points  on  the  left  and  right  sides  of  the  domain,  and  (5)  slabwise 
calculation  through  the  model  grid  to  improve  the  use  of  high-speed  cache  memory. 

The  sections  that  follow  include  (2)  a  description  of  the  model  physics,  (3)  a  description  of  the 
model  numerics,  (4)  a  discussion  of  the  limitations  of  the  model,  (5)  a  discussion  of  some  plans  for 
future  development  of  NCOM,  and  (6)  a  summary. 

2.  MODEL  PHYSICS 


2.1  Basic  Equations 


The  ocean  model  is  free  surface  and  employs  the  hydrostatic,  Boussinesq,  and  incompressible 
approximations.  The  equations  that  are  solved,  written  in  Cartesian  coordinates,  are 


T-r  =  -V  •  (vu)  +  Qu  +  fv  -  +  Fu  +  (km^)  , 


dt 


p0  dx 


dz 


(1) 


dv 

m 


1  dp  d  (  dv  \ 


(2) 


dp 

Tz  =  ~p9' 


du  dv  dw  _ 

vr  +  7T  +  IT  = 
dx  dy  dz 


f  -  -V  ■  (vT)  +QT  +  V,(A„VhT)  +  !  )  +  e-S’ 


dS 

dt 


^  q  /  Q  Qi  v 

\  =  -V  •  (vS)  +QS  +  Vh(AHVhS)  +  Q-z  (kh  —  j  , 


(3) 

(4) 

(5) 

(6) 
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p  =  p(T,S,z),  (7) 

where  t  is  the  time;  x.  y,  and  z  are  the  three  coordinate  directions;  u,  v,  and  w  are  the  three 
components  of  the  velocity  vector;  Q  is  a  volume  flux  source  term;  v  is  the  vector  velocity;  T 
is  the  potential  temperature;  S  is  the  salinity;  V/(  is  the  horizontal  gradient  operator;  /  is  the 
Coriolis  parameter;  p  is  the  pressure;  p  is  the  water  density;  p0  is  a  reference  water  density;  g 
is  the  acceleration  of  gravity;  Fu  and  Fv  are  horizontal  mixing  terms  for  momentum;  AH  is  the 
horizontal  mixing  coefficient  for  scalar  fields  (temperature  and  salinity);  Km  and  Kh  are  vertical 
eddy  coefficients  for  momentum  and  scalar  fields,  respectively;  Qr  is  the  solar  radiation;  and  7  is  a 
function  describing  the  solar  extinction. 

Density  p  must  be  calculated  from  T  and  S  using  a  suitable  equation  of  state.  Two  equations  of 
state  are  provided  within  NCOM,  the  polynomial  equation  of  state  of  Friedrich  and  Levitus  (1972) 
and  the  United  Nations  Educational,  Scientific,  and  Cultural  Organization  (UNESCO)  formula  as 
adapted  by  Mellor  (1991)  for  use  in  POM. 

2.2  Surface  and  Bottom  Boundary  Conditions 

Equations  (1)  to  (6)  are  subject  to  boundary  conditions  in  the  form  of  fluxes  and  stresses  at 
the  ocean’s  surface  and  bottom.  The  boundary  conditions  at  the  surface  at  z  —  £,  which  are  due 
to  fluxes  between  the  ocean  and  the  atmosphere,  are 


du_T * 

OZ  p0 

(8) 

dv  TV 

r~  5 

OZ  po 

(9) 

ts  &T  __  Qb  +  Qe  +  Qs 
oz  p0cp 

(10) 

f)Q 

KH-^  =  S\Z=C(EV-Pr), 

(11) 

where  tx  and  ry  are  the  x  and  y  components  of  the  surface  wind  stress,  Qb,  Qe,  and  Qs  are  the 
net  longwave  and  latent  and  sensible  surface  heat  fluxes,  Ev  and  Pr  are  the  surface  evaporation 
and  precipitation  rates,  and  cp  is  the  specific  heat  of  seawater.  At  the  ocean  bottom  at  2  =  FI.  the 
boundary  conditions  are 

KMT^=Cbu\v\,  (12) 


Km^z  =  CbV^’ 


(13) 


dT 

K^  =  0’ 


(14) 
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Ka%-  =  0.  (15) 

dz 

The  bottom  stress  is  parameterized  here  by  a  quadratic  drag  law  with  drag  coefficient  <‘i, ■  The 
value  of  Cfc  can  be  specified  or  calculated  in  terms  of  the  bottom  layer  thickness  A  Zb  and  the 
bottom  roughness  za  as 

q,  =  max 

where  n  =  0.4  is  von  Karman’s  constant,  and  is  a  minimum  value  for  q,.  This  expression  for 
Cfj  is  derived  by  assuming  a  logarithmic  boundary  layer  velocity  profile  near  the  bottom. 

2.3  Horizontal  Pressure  Gradient 


•  CK 


(16) 


Integration  of  the  hydrostatic  Eq.  (3)  in  the  vertical  from  a  depth  z  to  the  surface  at  z  =  ( 
yields 


f<  dp 

Jz  dz 


dz  =  -g 


The  vertical  pressure  gradient  integrates  exactly  to  give 


P( C)  -P(z)  -  ~9 


(17) 


(18) 


Taking  a  horizontal  derivative  (e.g.,  in  x)  and  using  Leibnitz’s  rule  for  the  differentiation  of  an 
integral, 


dp{ C)  _  dp 
dx  dx 


(19) 


Taking  p{()  ~  pn,  dividing  by  p0,  and  rearranging  terms,  the  expression  for  the  Boussinesq  hori¬ 
zontal  pressure  gradient  in  x  is 

L%.  =  LM1+S  +  Z  Op*.  (20) 

p0  OX  Po  OX  OX  p0  Jz  OX 


The  first  term  on  the  right  side  of  Eq.  (20)  is  the  pressure  gradient  at  the  surface  (i.e.,  the  surface 
atmospheric  pressure  gradient),  the  second  term  is  the  horizontal  pressure  gradient  due  to  differ¬ 
ences  in  the  surface  elevation,  and  the  third  term  is  the  horizontal  pressure  gradient  due  to  the 
density  field,  which  is  sometimes  referred  to  as  the  baroclinic  pressure  gradient.  The  horizontal 
pressure  gradient  in  y  is  calculated  in  similar  fashion. 


2.4  Horizontal  Mixing 

The  model  uses  the  Laplacian  form  of  horizontal  mixing,  and  two  horizontal  mixing  parame- 
terizations  are  provided.  One  is  based  on  maintaining  a  maximum  horizontal  grid-cell  Reynolds 
number,  and  the  other  is  the  Smagorinsky  scheme  (Smagorinsky  1963),  which  is  used  in  POM. 

For  the  grid-cell  Reynolds  number  scheme,  the  horizontal  friction  terms  in  the  momentum 
equations  use  the  form 
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d(  dv\  d  (  .  dv 
dx  \M~dx)  +  fry  \AMdy 


(22) 


The  horizontal  mixing  coefficient  Am  is  set  equal  to  the  maximum  of  a  constant  background  value 
A0  and  a  value  needed  to  keep  the  grid-cell  Reynolds  Number  Reg  below  a  maximum  specified 
value,  i.e.,  in  the  ^-direction, 


Am  =  max 


|ti|Aa: 


R, 


9  J 


(23) 


and  in  the  //-direction 


Am  =  max 


(24) 


where  Aar  and  Ay  are  the  horizontal  grid  spacing  in  x  and  y.  respectively.  The  value  of  Rtg  is 
typically  set  to  a  value  in  the  range  of  10  to  100. 


For  the  Smagorinsky  mixing  scheme,  the  horizontal  friction  terms  have  the  form 


The  horizontal  mixing  coefficient  Am  is  calculated  as  a  function  of  the  local  horizontal  grid  resolu¬ 
tion  and  velocity  shear,  i.e., 

<  „  \  a  ( fdu\2  1  ( dv  ®u\ 2  fdv\2\2  .  , 

Am  =  CsmagAxAy  +  -  (—  +  +  [^)  J  ■  (  7) 

The  magnitude  of  the  Smagorinsky  eddy  coefficient  calculation  is  scaled  by  the  constant  Csmag. 
Values  used  for  Csmag  range  from  about  0.02  to  0.5.  Large  values  of  Csmag  tend  to  dissipate  smaller 
features,  whereas  values  that  are  too  small  can  lead  to  excessive  numerical  noise  and/or  instability. 
A  typical  value  that  is  used  is  0.1. 


The  grid-cell  Reynolds  number  scheme  is  simpler  and,  with  the  mixing  coefficients  scaled  ac¬ 
cording  to  the  velocity  magnitude,  is  specifically  geared  toward  suppressing  noise  generated  by 
numerical  advection.  The  Smagorinsky  scheme  scales  the  rate  of  mixing  according  to  the  hori¬ 
zontal  velocity  shear  and  might  be  considered  to  be  more  physically  based.  The  eddy  coefficients 
calculated  by  the  Smagorinsky  scheme  are  isotropic  and  are  independent  of  coordinate  rotation; 
those  calculated  by  the  grid-cell  Reynolds  number  method  are  not. 

The  form  used  for  horizontal  diffusion  of  T  and  S  is  the  Laplacian  form  expressed  in  Eqs.  (5) 
and  (6).  NCOM  provides  for  allowing  Ah  to  be  some  fraction  of  Am  via  specification  of  an  inverse 
Prandtl  Number.  Setting  A //  =  Am  is  a  common  choice,  but  for  some  applications  it  may  be 
desirable  to  set  Ah  different  from  Am- 
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2.5  Vertical  Mixing 


The  vertical  eddy  coefficients  are  specified  as 


Km  =  max 


Km0.  A' a/,  ,  KM2 , 


w\Az' 

~r^t\  5 


(28) 


Kji  =  max 


Kho-  KHi,Kh2, 


w\Az 

Re. 


(29) 


where  Km0  and  Kh0  are  small,  background  values,  K a/,  and  Ku1  are  calculated  from  a  turbu¬ 
lence  model,  Km0  and  Kh2  parameterize  turbulent  mixing  by  unresolved  processes  at  near-critical 
Richardson  Numbers  (Large  et  al.  1994),  Az  is  the  vertical  grid  spacing,  and  Re,  is  a  maximum 
grid-cell  Reynolds  number  for  the  vertical  direction. 


The  background  values  K\!o  and  Kf[()  are  constants.  Their  purpose  is  to  parameterize  weak, 
vertical  mixing  processes  that  are  not  accounted  for  by  the  other  mixing  parameterizations.  They 
are  generally  kept  quite  small  since  large  values  would  be  unrealistic  and  would  excessively  erode 
the  stratification  and  damp  the  flow,  especially  in  shallow  water.  Typical  values  are  10  5  m2/s  or 
less. 


Two  turbulence  schemes  are  provided  for  calculating  K and  Kj{x ,  the  MYL2.5  scheme  (Mellor 
and  Yamada  1982;  Mellor  1996),  which  is  used  in  POM,  and  the  simpler  MYL2  scheme  (Mellor  and 
Yamada  1974;  Mellor  and  Durbin  1975).  In  both  of  these  schemes,  Kmx  and  Kh1  are  calculated  as 

KMl  =  tqSM ,  (30) 


KHl  =  eqSH ,  (31) 

where  i  is  a  vertical  turbulent  length  scale,  \ql  is  the  turbulent  kinetic  energy  (TKE),  and  S\i  and 
SH  are  stratification  functions  that  describe  the  effect  of  stratification  on  the  vertical  mixing. 


The  MYL2.5  scheme  provides  a  prognostic  equation  to  calculate  the  TKE  that  includes  advec- 
tive  and  diffusive  transport.  Another  prognostic  equation  for  the  quantity  q2f  is  used  to  provide 
an  estimate  for  the  vertical  turbulence  length  scale  i.  These  two  equations  are 


dq 2 
dt 


-V  •  (vq2)  +  Qq 2  +  V/l(Af/V/lg2)  +  ^ 


dq 2 
dz 


+  2Kj\[1 


+  2  ^-KHl 

Po 


2  q3 

M’ 


dq2d 

~df 


-V  •  (vq2t)  +  QQ2?  +  Vh{AHVhq2l)  + 


+  E 


(32) 


(33) 


where 
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W  =  1  +  E2 


(34) 


L  1  =  ((  —  z  +  zs)  l  +  (z-H  +  z0)  \ 


(35) 


dp  dp  n 
d^  =  d^^Cs  ~dz 


(36) 


In  the  above  equations,  Kqi  is  the  vertical  diffusion  coefficient  for  q2  and  q2l.  which  is  taken  to  be 
proportional  to  KHl  (Kqi  =  0A1KHi),  zs  is  the  surface  roughness,  cs  is  the  speed  of  sound,  and 
&i,  Ei,  E-2,  and  E$  are  constants  (see  Table  2). 


As  discussed  by  Mellor  and  Yamada  (1982),  Eq.  (33),  which  is  used  to  obtain  the  turbulence 
length  scale,  is  somewhat  ad  hoc.  Most  of  the  methods  used  to  estimate  turbulence  length  scales 
resort  to  some  degree  of  empiricism. 

Mellor  and  Yamada  (1982)  refer  to  W  as  a  “wall-proximity”  function,  which  is  used  to  scale 
l  near  the  surface  and  bottom.  This  function  has  been  modified  from  the  form  used  in  POM  to 
account  for  the  surface  roughness  length  zs  and  the  bottom  roughness  length  z0  (in  POM,  L  in 
Eq.  (35)  tends  to  zero  at  the  top  and  bottom).  Since  zs  and  z0  are  usually  fairly  small  relative  to 
the  vertical  grid  resolution,  the  effect  of  this  change  is  not  usually  significant.  However,  in  the  case 
of  strong  winds  and  breaking  waves  at  the  surface,  the  surface  roughness  can  significantly  increase 
the  mixing  in  the  surface  mixed  layer  (Craig  and  Banner  1994;  Craig  1996). 

Table  2  —  Constants  for  MYL2.5  Turbulence  Model 


parameter 

value 

ai 

0.92 

a,2 

0.74 

h 

16.6 

b2 

10.1 

Cl 

0.08 

Ei 

1.8 

i?2 

1.33 

1.0 

The  density  p,  which  is  used  to  calculate  the  vertical  buoyancy  gradient  in  the  TKE  equation, 
must  not  include  the  effect  of  local  changes  in  pressure  on  the  density;  otherwise,  the  vertical 
stability  will  be  overestimated.  If  p  is  calculated  as  an  in  situ  density,  where  the  effect  of  pressure 
on  the  density  is  accounted  for,  Eq.  (36)  provides  the  correction  to  remove  the  effect  of  local  pressure 
changes  on  the  vertical  density  gradient.  If  the  model  density  does  not  include  the  effect  of  pressure, 
one  can  set  p  =  p  (in  shallow  water,  the  effects  of  pressure  on  the  density  are  small  and  are  often 
neglected). 

Prognostic  Eqs.  (32)  and  (33)  require  boundary  conditions  at  the  surface  and  bottom.  At  the 
surface  at  z  =  (,  the  boundary  conditions  are  given  by 
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NCOM  also  provides  the  option  of  specifying  the  flux  of  TKE  at  the  surface  rather  than  the 
value  itself.  The  surface  TKE  flux  is  currently  parameterized  in  terms  of  the  surface  wind  stress, 
i.e.,  at  z  —  ( 

KM^£  =  cwave^y,  (39) 

where  cwave  is  a  constant  that  scales  the  TKE  input  from  the  waves.  (Note  that  if  wave  data  are 
available,  it  might  be  preferable  to  parameterize  the  surface  TKE  flux  in  terms  of  the  actual  sea 
state.)  This  boundary  condition  for  the  surface  TKE  flux,  along  with  accounting  for  the  surface 
roughness  due  to  surface  waves,  allows  for  the  simulation  of  the  wave-mixing  layer  near  the  surface 
in  which  there  is  enhanced  mixing  from  the  energy  input  from  breaking  waves  (Craig  and  Banner 
1994). 

At  the  bottom  at  z  =  H,  the  boundary  conditions  are 

q2  =  bjcb{u2  +  w2)V  (40) 

q2i  =  bfcbiu2  +  v2)ikz0.  (41) 


The  stratification  functions  Sjq  and  Sh  are  calculated  as 


Sh 


Ci 

1  -C2Gh' 


(42) 


_  C3  +  C4SH 
~  1  -  C,Gh  ’ 


(43) 


where 


Gh  —  min 


0.028, 


l2g  dp 


(44) 


q2Po  dz 

The  constants  C\  -  C-,  are  calculated  from  the  basic  turbulence  constants  (ai,  a2,  b\,  b2,  and  c\)  as 


C\  —  ^2(^1  “  6ai)/6x, 

(45) 

C2  —  ( 1 801  +  3^2)? 

(46) 

C3  =  a1(61(l-3C1)-6a1)/&i, 

(47) 
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C4  —  ai(18ai  +  9^2), 


(48) 


C5  =  9aia2.  (49) 

Table  2  lists  the  values  of  the  constants  used  in  the  MYL2.5  turbulence  model. 


The  MYL2  turbulence  model  assumes  that  there  is  approximately  a  local  balance  between 
shear  production,  buoyancy  production,  and  dissipation  in  the  TKE  equation  (these  are  the  last 
three  terms  in  Eq.  (32)).  With  this  assumption,  q  can  be  calculated  algebraically  from  the  mean 
vertical  density  and  velocity  gradients 


<?  = 


+  s„±f 

Po  OZ. 


1 

2 


(50) 


An  empirical  formulation  is  used  to  estimate  £.  The  turbulent  length  scale  l  is  set  to  zero 
outside  of  turbulent  layers,  and  within  a  turbulent  layer  is  scaled  as  the  vertical  distance  to  the 
edge  of  the  turbulent  layer,  i.e.,  if  zt  is  the  distance  to  the  top  of  the  turbulent  layer  and  zt>  is  the 
distance  to  the  bottom  of  the  turbulent  layer,  then  the  local  turbulent  length  scale  £  within  the 
turbulent  layer  is  calculated  as 


i  =  1  +  Zb1)  \  (51) 

In  the  surface  mixed  layer,  the  surface  roughness  is  added  to  zt,  and  in  the  bottom  boundary  layer 
the  bottom  roughness  is  added  to  zb  so  that  the  value  of  i  reflects  the  surface  and  bottom  roughness. 
Calculated  from  Eq.  (51),  i  has  a  quadratic  profile  within  a  turbulent  layer,  and  the  value  of  £  at 
the  center  of  a  turbulent  layer  is  about  10%  of  the  thickness  of  the  layer  (for  small  values  of  zs  and 
zQ ).  This  is  a  somewhat  crude  parameterization  of  £,  but  this  is  mitigated  somewhat  by  the  fact 
that  the  depth  of  turbulent  layers  predicted  by  the  MYL2  scheme  in  density  stratified  conditions 
tends  to  be  only  weakly  dependent  on  the  value  of  t.  There  is  also  the  point  that  more  elaborate 
turbulence  parameterizations  frequently  do  not  describe  the  details  of  turbulent  mixing  very  well 
since  they  do  not  account  for  or  adequately  resolve  some  of  the  significant  processes. 

The  stratification  functions  Sm  and  Sh  for  MYL2  are  calculated  from  the  Richardson  Number 
Ri  as 


SH  =  C'1-  C'2R, 


(52) 


Sm  =  Sh 


c3  -  c'ar 

C's-QR' 


where 


R  = 


Rf 

1  -Rf' 


(53) 


(54) 


Rf  =  C '(•  +  C$Rj  —  ((Cg  Ri  +  C 9 ) /?,<  +  C72) 2 , 


(55) 


Description  of  NCOM  Version  1.0 


1  11 

C'l  =  a2(b  1  -  6ai)/6i, 

(56) 

C'2  =  02(1801  +  3&2)/&i, 

(57) 

C3  =  ai(6i(l  -  3ci)  -  6oi)/o2, 

(58) 

C\  =  oi(18oi  +  9a2)/a2, 

(59) 

C5  =  b\  —  601, 

(60) 

Cq  —  9ai  +  352, 

(61) 

1  C*3 

7  2  c^  +  cy 

(62) 

r,  1  CJ  +  Ce 

8  2  C'3  +  C\  ’ 

(63) 

nt 

S~i!  0  ^5 

O9  —  C3  +  C4  # 

(64) 

In  these  equations,  Rf  is  a  flux  Richardson  Number  (Rj  —  R,Kn /Km  =  R.,Sh  /  Sm)-,  and  cq,  0,2,  h, 
b2,  and  ci  are  the  the  same  basic  turbulence  constants  as  used  for  the  MYL2.5  model.  Table  3  lists 
the  values  for  the  constants  used  for  the  MYL2  turbulence  model.  Note  that  the  basic  turbulence 
constants  used  here  for  the  MYL2  model  are  as  used  by  Mellor  and  Durbin  (1975)  and  are  slightly 
different  from  the  values  used  in  the  MYL2.5  model. 

Table  3  —  Constants  for  MYL2  Turbulence  Model 


parameter 

value 

a\ 

0.78 

a2 

0.78 

bi 

15.0 

b2 

8.0 

c  1 

0.056 

The  MYL2  turbulence  model  as  described  here  was  compared  with  the  MYL2.5  scheme  as 
used  in  POM  for  some  simple  cases  of  surface-layer  mixing  by  winds  and  heat  fluxes  and  bottom- 
layer  mixing  by  tidal  currents  and  was  found  to  give  similar  turbulence  mixing  scales  and  similar 
turbulent  layer  depths  for  these  cases  (Martin  1985;  Martin  1986;  Martin  et.  al.  1998).  Because  of  its 
simplicity,  the  MYL2  turbulence  model  is  relatively  efficient  and  does  not  carry  the  computational 
burden  of  the  MYL2.5  scheme,  which  requires  the  model  to  carry  and  calculate  two  additional 
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prognostic  variables.  However,  the  MYL2.5  scheme  has  some  advantages.  For  example,  in  high- 
resolution  simulations,  transport  of  TKE  will  be  more  significant  and  the  assumption  of  local 
equilibrium  of  turbulence  by  the  MYL2  model  may  be  less  accurate.  Also,  because  it  does  not 
account  for  vertical  diffusion  of  TKE,  the  MYL2  model  cannot  simulate  the  wave  mixing  layer  that 
occurs  near  the  surface  in  strong  winds. 

Tests  of  upper-ocean  mixing  in  the  open  ocean  with  turbulence  models  such  as  MYL2  and 
MYL2.5  have  frequently  found  that  the  models  do  not  predict  as  much  mixing  as  is  observed 
(Martin  1985;  Large  et  al.  1994;  Kantha  and  Clayson  1994).  It  has  been  hypothesized  that  this 
is  due,  not  so  much  to  the  fact  that  the  mixing  mechanisms  in  these  models  are  incorrect,  but 
that  many  mixing  processes  are  not  accounted  for  in  these  simulations,  including  many  sources  of 
background  shear  (internal  waves,  inertial  gravity  wave  pumping,  etc.)  and  Langmuir  circulations. 

NCOM  provides  an  option  to  include  the  vertical  mixing  enhancement  scheme  of  Large  et 
al.  (1994).  The  purpose  of  this  mixing  scheme  is  to  account  for  unresolved  mixing  processes  by 
extending  the  mixing  of  typical  oceanic  turbulence  models  above  the  normal  critical  Richardson 
number  value  of  0.2  to  0.5.  The  Large  et  al.  enhancement  scheme  extends  the  mixing  to  Ri  =  0.7 
and  is  described  by 

K0  Ri  <0 

KM2=Kh2=  K0(l  -  (Ri/0.7)2)3  0  <Ri  <0.7 

0  Ri  >  0.7, 

where  K0  =  50  cm2/s.  This  scheme  was  utilized  by  Large  et  al.  in  conjunction  with  an  adaptation 
of  the  atmospheric  boundary  layer  model  of  Troen  and  Mahrt  (1986)  and  by  Kantha  and  Clayson 
(1994)  in  conjunction  with  the  MYL2.5  turbulence  closure  model.  Both  Large  et  al.  and  Kantha 
and  Clayson  found  that  the  addition  of  this  mixing  improved  agreement  of  predictions  of  the  ocean 
surface  mixed  layer  with  observations  from  several  open-ocean  data  sets.  However,  it  is  not  clear 
whether  such  a  mixing  enhancement  is  needed  in  shallow,  coastal  water. 

2.6  Vertically  Averaged  Equations 

The  depth-averaged  momentum  and  continuity  equations  are  needed  to  calculate  the  free- 
surface  mode  in  the  model.  Substituting  the  form  of  the  horizontal  pressure  gradient  from  Eq.  (20) 
into  the  momentum  Eqs.  (1)  and  (2)  and  integrating  from  the  bottom  to  the  surface  gives  the  form 
of  the  momentum  equations  used  for  the  calculation  of  the  free-surface  mode 


dDv 

dt 


(66) 


where  D  is  the  total  depth  (D  =  £  —  H),u  and  v  are  the  depth-averaged  horizontal  velocities,  and 
Q  is  the  depth-averaged  mass  flux  source  term,  and 


G 


u  — 


— V  •  (vu)  +  Qu  +  fv  — 


1  dpjQ 

Po  dx 


9_ 

Po 


l 


£  dp  ,  d 

~^~dz  +  Fu  +  — 
ox  oz 


(67) 
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Gv  =  -V  •  (vu)  +  Qv  —  fu  - 


1  dp{  0 
po  dy 


9_ 

Po 


L 


^dz  +  Fv 

dy 


(68) 


The  depth-integrated  continuity  equation  is 


dC  d(Du)  d[Dv)  J, 

at  dx  ay  w' 


3.  MODEL  NUMERICS 
3.1  Horizontal  Grid 

The  model  uses  an  orthogonal,  curvilinear,  horizontal  grid  as  used  in  POM,  which  allows 
adaptation  of  the  grid  to  different  map  projections  and  also  allows  the  grid  to  be  set  up  with  some 
mild  curvature.  The  form  of  the  model’s  equations  in  orthogonal,  curvilinear,  horizontal  coordinates 
is  presented  in  Appendix  A.  As  a  practical  matter,  the  main  differences  from  the  Cartesian  equations 
with  constant  grid  spacing  in  x  and  y  with  regard  to  converting  the  equations  to  finite  difference 
form  are  (1)  the  horizontal  grid-cell  dimensions  in  x  and  y.  Ax  and  Ay,  respectively,  can  vary 
spatially  and  must  be  stored  as  two-dimensional  (2-D)  arrays,  (2)  the  fluxes  between  grid  cells  must 
account  for  the  changing  size  of  the  grid  cells,  and  (3)  correction  terms  are  needed  to  account  for  the 
exchange  between  u  and  v  momentum  due  to  horizontal  transport  along  curving  grid  coordinates. 

The  spatial  finite  differencing  in  NCOM  is  done  in  conservative  form  with  the  advective  and 
diffusive  transport  between  grid  cells  calculated  as  fluxes  between  the  grid  cells.  The  use  of  this 
form  maintains  conservation  of  the  scalar  model  fields  for  transport  between  grid  cells  regardless 
of  how  the  size  of  the  grid  cells  change. 

With  transport  conservation  satisfied  by  the  flux  form  used  for  differencing,  the  main  remaining 
requirement  for  the  use  of  a  curvilinear  grid  is  that,  since  the  horizontal  velocity  components  u  and 
v  follow  (i.e.,  are  directed  along)  curving  horizontal  coordinates,  there  is,  in  effect,  a  conversion 
of  momentum  between  the  two  horizontal  velocity  components  as  momentum  is  transported  along 
a  curving  grid  coordinate.  There  is  a  correction  for  both  the  horizontal  advection  and  diffusion 
terms  to  account  for  interchange  between  u  and  v  momentum  on  a  curvilinear  grid.  POM  does  not 
provide  a  curvature  correction  for  the  horizontal  diffusion  terms  in  the  momentum  equations  and 
none  is  currently  provided  in  NCOM.  The  assumption  is  that  the  transport  of  momentum  due  to 
horizontal  diffusion  is  sufficiently  small  that  error  due  to  neglect  of  the  curvature  correction  to  this 
term  will  not  be  significant. 

To  keep  truncation  errors  associated  with  curvilinear  grids  relatively  small,  a  rule  of  thumb  in 
using  such  grids  is  to  not  change  the  size  of  the  grid  spacing  by  more  than  about  10%  between  suc¬ 
cessive  grid  cells.  With  the  simple  two-point  averages  and  differences  used  for  the  finite  differencing, 
the  accuracy  of  spatial  interpolations  and  gradients  becomes  first-order  rather  than  second-order 
if  the  change  in  size  between  successive  grid  cells  becomes  more  than  a  small  fraction  of  the  grid 
spacing. 

The  philosophy  that  has  been  followed  to  date  in  applying  NCOM  has  been  to  avoid  using 
grid  curvature,  except  as  needed  to  adapt  the  model  grid  to  large-scale  map  projections  where  the 
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grid  curvature  is  very  gradual,  to  minimize  spatial  truncation  errors.  In  our  coastal  simulations  to 
date,  we  have  not  encountered  a  need  to  use  strong  grid  curvature.  The  entire  grid  can  be  rotated 
to  provide  a  desired  orientation  along  a  section  of  coastline,  and  shrinkwrapping  can  be  used  to 
eliminate  calculations  over  land  areas  along  the  sides  of  the  domain.  Most  coastlines  are  so  irregular 
that  curvilinear  coordinates  cannot  be  curved  sufficiently  to  follow  them  with  any  fidelity. 

Figure  1  shows  the  horizontal  arrangement  of  some  of  the  variables  on  the  model  grid.  The 
horizontal  arrangement  of  the  variables  uses  the  form  of  the  Arakawa,  C  grid.  With  the  C  grid, 
scalar  fields  such  as  T,  S ,  and  p  are  located  at  the  grid-cell  midpoints,  and  each  of  the  velocity 
components  is  located  at  the  center  of  the  grid-cell  face  to  which  it  is  normal. 


Fig.  1  —  Horizontal  layout  of  variables  on  model  grid 


3.2  Vertical  Grid 

The  model  uses  a  combined  sigma/;?- level  vertical  grid  with  sigma  layers  near  the  surface  and 
^-levels  below  a  depth  that  can  be  specified  by  the  user. 

Figure  2  illustrates  the  different  ways  the  combined  sigma/ 2-level  vertical  grid  can  be  set  up. 
Figure  2(a)  shows  the  vertical  grid  set  up  with  a  single  sigma  layer  at  the  surface  and  with  2-levels 
below.  Since  the  model  has  a  free  surface,  at  least  one  sigma  layer  is  needed  at  the  surface  to  allow 
for  changes  in  the  surface  elevation. 

If  the  changes  in  the  surface  elevation  are  large  relative  to  the  grid  resolution  desired  near  the 
surface,  a  single  sigma  layer  may  not  be  sufficient  to  resolve  the  changes  in  the  surface  elevation.  In 
this  case,  several  sigma  layers  can  used,  and  the  changes  in  the  surface  elevation  will  be  distributed 
among  them  (Fig.  2(b)). 

If  the  water  depth  becomes  shallower  than  the  depth  that  defines  the  transition  from  sigma 
layers  to  2-levels  (2ff),  the  sigma  layers  will  shallow  uniformly  as  the  bottom  depth  decreases.  Figure 
2(c)  shows  a  grid  in  which  the  sigma  layers  extend  to  the  bottom  in  the  shallow  water  on  the  shelf, 
and  2-levels  are  used  in  the  deeper  water  off  the  shelf.  Figure  2(d)  shows  a  grid  in  which  sigma 
layers  are  used  all  the  way  to  the  bottom  everywhere. 
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Fig.  2  —  Illustration  of  the  different  ways  the  combined  sigma/2-level  grid  can  be  set  up:  (a)  with  a  single  sigma  layer 
at  the  surface,  (b)  with  several  sigma  layers  at  the  surface  to  resolve  surface  elevation  changes,  (c)  with  sigma  layers 
in  the  shallow  water  and  2-levels  in  the  deep  water,  and  (d)  with  sigma  layers  all  the  way  to  the  bottom  everywhere 


Figure  3  shows  the  vertical  arrangement  of  some  of  the  model  variables  on  the  model  grid.  As 
in  the  layout  of  the  horizontal  grid,  the  main  scalar  fields  are  located  at  grid-cell  centers,  and  the 
velocity  components  are  located  at  the  center  of  the  grid-cell  face  to  which  they  are  normal. 


The  coordinate  transformation  for  the  sigma  coordinate  part  of  the  grid  is  given  by 

*-< 


a  = 


Da 


(70) 


where  Da  —  Q  -  ma x{H,za).  Hence,  a  varies  from  a  —  0  at  the  free  surface  to  ct  =  —  1  at  the 
bottom  interface  of  the  lowest  sigma  layer.  Each  sigma  layer  is  a  fixed  fraction  of  the  total  depth 
of  the  sigma  grid  Da . 


Similar  to  the  implementation  of  orthogonal,  curvilinear  coordinates,  the  implementation  of 
sigma  coordinates  is  primarily  a  matter  of  accounting  for  the  changing  vertical  thickness  of  the 
layers  in  calculating  fluxes  between  adjacent  grid  cells  and  in  calculating  changes  within  the  grid 
cells.  Note  that  with  sigma  coordinates,  the  changes  in  the  thickness  of  the  grid  cells  occurs  not 
only  horizontally  within  a  layer  but  also  from  timestep  to  timestep  because  of  the  changing  surface 
elevation.  To  avoid  having  to  use  a  large  number  of  3-D  arrays  to  store  the  thickness  of  the  grid 
cells  at  different  locations  and  time  levels,  the  grid  thickness  Az  for  a  sigma  layer  is  expressed  as 
the  product  of  the  fractional  thickness  of  the  sigma  layer  Act  times  the  depth  from  the  surface  to 
the  bottom  of  the  sigma  coordinate  grid,  i.e.,  as 


Az  =  A  <jDa. 


(71) 
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Fig.  3  —  Vertical  layout  of  variables  on  model  grid 


The  vertical  velocity  on  the  sigma  coordinate  grid  uj  is  the  vertical  velocity  relative  to  the  sigma 
surfaces.  This  velocity  is  not  the  true  vertical  velocity  since  it  is  missing  the  vertical  component 
of  the  flow  along  the  sigma  layers  (if  they  are  sloping)  as  well  as  the  vertical  velocity  due  to  the 
vertical  motion  of  the  sigma  surfaces  themselves  caused  by  the  rise  and  fall  of  the  sea  surface.  The 
velocity  u>  is  related  to  the  Cartesian  vertical  velocity  w  as 

w  =  w  +  (1  -  <t)%  +  u#(C  +  oDa)  +  v^-(C  +  (72) 

dt  dx  oy 

Since  the  surfaces  on  the  z-level  grid  are  level  and  fixed  in  time,  w  ~  u;  on  the  £-level  grid. 


The  form  of  the  model  equations  in  sigma  coordinates  is  presented  in  Appendix  B.  The  only 
significant  modification  to  the  basic  equations  for  the  changing  depth  of  the  sigma  layers,  with 
regard  to  converting  the  equations  to  finite  difference  form,  is  a  correction  for  the  horizontal  pressure 
gradient.  The  horizontal  pressure  gradient  calculation  in  sigma  coordinates  contains  an  extra  term 
to  remove  the  vertical  change  in  pressure  between  neighboring  points  within  a  sigma  layer,  so  that 
the  net  pressure  change  that  is  computed  will  be  approximately  along  a  level  surface  (Blumberg 
and  Mellor  1987).  The  form  of  the  horizontal  pressure  gradient  in  sigma  coordinates  is  modified 
from  the  Cartesian  form  (Eq.  (20))  used  on  the  z-level  part  of  the  grid  to 


1  dp  _  1  dp{ C)  3C  ,  f°  \dp< 
pQdx  p0  dx  ^  dx  p0  J (j  dx° 


o  dDa  dp ' 
Da  dx  do  _ 


do , 


(73) 


where  the  term  ^\a  is  the  density  gradient  taken  along  a  surface  of  constant  o. 


A  potential  problem  with  this  calculation  of  the  pressure  gradient  in  sigma  coordinates  is  that 
the  vertical  component  of  the  pressure  change  along  a  sloping  sigma  surface  is  frequently  much 
larger  than  the  horizontal  component  along  a  level  surface  (Haney  1991).  In  this  case,  the  desired 
horizontal  component  is  calculated  as  the  small  difference  between  two  large  terms  and  is  subject  to 
significant  truncation  error.  An  expedient  that  is  commonly  used  to  reduce  this  error  is  to  subtract 
the  horizontally  averaged  density  profile  from  the  3-D  density  field  when  calculating  the  horizontal 
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pressure  gradient,  so  that  the  main  component  of  the  vertical  change  in  pressure  is  removed  from 
the  calculation,  i.e.,  p  in  Eq.  (73)  is  replaced  by  p  —  p(z),  where  p(z)  is  the  horizontal  mean  density 
profile  (Blumberg  and  Mellor  1987). 

Strictly  speaking,  for  a  full  transformation  of  the  equations  to  sigma  coordinates,  the  horizontal 
diffusion  terms  should  be  corrected  for  the  transformation,  so  that  they  will  still  represent  diffusion 
along  level  surfaces.  However,  NCOM  (and  POM)  use  the  approximation  discussed  by  Mellor 
and  Blumberg  (1985),  who  argued  that  diffusion  along  the  sigma  surfaces  rather  than  along  level 
surfaces  was,  in  general,  more  appropriate  for  sigma  coordinate  models,  particularly  for  proper 
simulation  of  the  bottom  boundary  layer. 

However,  in  regions  where  there  are  large  changes  in  the  bottom  depth,  the  horizontal  diffusion 
along  sloping  sigma  layers  can  cause  severe  cross- isopycnal  diffusion  (Paul  1994).  An  expedient  that 
is  sometimes  used  with  sigma  coordinate  models  (Mellor  and  Blumberg  1985),  and  is  an  option  in 
NCOM,  is  to  subtract  a  smooth  background  field  from  the  T  or  S  fields  when  calculating  horizontal 
diffusion.  By  calculating  the  horizontal  diffusion  based  on  the  anomaly  from  a  smooth  background 
field,  most  of  the  component  of  vertical  diffusion  that  occurs  when  diffusion  is  calculated  along 
sloping  sigma  layers  is  eliminated. 

In  domains  where  the  T  and  S  fields  don’t  vary  much,  the  background  T  and  S  fields  can  be 
calculated  as  the  horizontally  averaged  T  and  S  profiles.  An  alternative  strategy  is  to  use  smooth 
but  horizontally  varying  background  fields  to  accommodate  changes  in  the  structure  of  the  T  and 
S  fields  in  different  parts  of  the  model  domain.  The  background  fields  can  also  be  periodically 
updated  to  accommodate  changes  in  T  and  5  that  occur  in  time.  The  use  of  these  procedures  can 
significantly  reduce  the  problem  of  severe  cross-isopycnal  diffusion  (however,  as  discussed  in  Section 
4.2,  they  can  introduce  other  problems). 

On  the  2-level  part  of  the  grid,  the  bathymetry  is  rounded  to  the  nearest  2-level.  This  is  the 
simplest  way  to  implement  bathymetry  in  a  2-level  model  and  is  the  way  the  bathymetry  has  been 
incorporated  into  a  number  of  ocean  models,  including  the  various  versions  of  the  Bryan-Cox  model 
(Bryan  1969;  Killworth  et  al.  1991;  Dukowicz  and  Smith  1994),  the  Haney  model  (Haney  1974), 
and  the  DieCAST  model  (Dietrich  and  Ko  1994).  There  are,  however,  a  number  of  limitations  to 
this  representation  of  the  bathymetry,  which  are  discussed  in  Section  4.3. 

3.3  Spatial  Differencing 


Spatial  interpolations  and  gradients  use  second-order  centered  averages  and  differences.  With 
second-order,  centered  interpolations,  the  value  of  a  variable  ^  at  a  location  between  points  at 
which  the  variable  is  defined  is  evaluated  as  the  average  of  the  values  on  either  side,  e.g.,  for  an 
interpolation  in  the  x  direction, 


<t>  ~  x(^a;+Aa;/2  +  0 x-Ax /?)■ 


With  second-order,  centered  differencing,  the  gradient  of  (j>  at  x  is  calculated  as  the  difference  of 
the  values  on  either  side 


d<j) 

Ox 


(75) 
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3.4  Temporal  Differencing 


The  leapfrog  scheme  is  used  for  temporal  differencing.  The  temporal  differencing  of  the  3-D 
equations  will  be  illustrated  here  with  just  the  u  momentum  equation  since  the  treatment  of  the 
other  model  variables  is  similar.  In  the  following  discussion,  n  will  denote  model  values  at  the 
current  time  level  (i.e.,  values  calculated  on  the  previous  timestep),  n  +  1  will  denote  the  newly 
calculated  values,  and  n  -  1  will  denote  values  at  the  previous  time  level. 


With  the  leapfrog  scheme,  most  of  the  terms  (i.e.,  the  advection,  baroclinic  pressure  gradient, 
and  Coriolis  terms)  are  centered  in  time  at  n.  The  horizontal  diffusion  terms  are  evaluated  at  the 
n  —  1  time  level  for  the  variable  being  diffused  (since  evaluation  of  the  variable  being  diffused  at  the 
central  time  level  of  a  leapfrog  scheme  is  numerically  unstable) ,  and  the  vertical  diffusion  terms  are 
treated  implicitly  so  as  to  avoid  the  timestep  restriction  for  explicit  vertical  diffusion  (the  high  rates 
of  vertical  diffusion  that  are  frequently  calculated  by  the  turbulence  schemes  would  require  a  very 
small  timestep  for  stability  if  the  vertical  diffusion  were  explicit).  Hence,  the  temporal  differencing 
of  the  u  momentum  equation  is  of  the  form 


fatu 

2A  t 


un+1  —  un  1 
2A t 


=  -V  •  (vu)n  +  fvr‘ 


1  dv11 

i  jpn- 
o  "b  xu 
Po  OX 


+ 


d_ 

dz 


dun+l \ 

dz  )  ' 


(76) 


An  Asselin  filter  is  used  to  suppress  the  time  splitting  that  can  occur  with  leapfrog  (Asselin 
1972).  The  Asselin  filter  is  applied  to  the  model  fields  at  time  level  n,  after  the  new  values  at  n  +  1 
have  been  calculated,  by  averaging  in  a  bit  of  the  values  at  the  n  +  1  and  n  —  1  time  levels,  i.e., 

4>n  =  v((l)n+1 +  <t>n-1)  +  (l-2u)(j)n,  (77) 

where  v  is  the  filter  coefficient.  If  Eq.  (77)  is  rewritten  as 

(f)n  =  (f/1  +  £/(</>n+1  —  2cj)n  +  (78) 

the  filter  looks  like  a  numerical  diffusion  term,  which  is  the  way  that  it  behaves.  A  typical  value 
used  for  v  is  0.05. 


3.5  Finite  Difference  Form  of  the  Model  Equations 


As  noted  earlier,  the  model  equations  are  finite  differenced  in  flux-conservative  form.  The  full 
finite  difference  form  of  the  basic  3-D  model  equations  used  in  NCOM  is 


AxuAyu 

2At 


1 


S2t(Azuu)  =  - AyuAzugSx(C  +  Catm  -  C tp)  ~  AyuAzu-6x(Pi) 


Po 


+  AxAyAz(f  +  Ccurv )& 

—  6x(AyuAzuuaXTix)  —  Sy(Axv  AzvvaXuy) 

—  5z(AxAyu>xuz)  +  AxAyAzQxusor  +  F* 


Km 3 


(Azw  )n+1 
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,71+ 1 


+  AxuA  yu8z 


(79) 
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A„V  An,V  1 

tA;7  S2t(Azvv)  =  -AxvA zvgSy(C  +  Cat m  ~  Ctp)  -  AxvAzv-6y(Pi ) 

2  At  Po 

—  AxAyAz(f  +  CCUrv)uxV 

—  8X  ( Ayu  AzuuaVvx )  —  8y(Axv  AzvvaXvy) 

—  8z(AxAyujyvz)  +  AxAyAzQIJvsor  +  F* 

+*x’av’5*(je^s*v"+')  m 

^^82t{Az)  =  -5x(AyuAzuua)  +  -8y(AxvAzvva)  +  -6z(AxAyuj) 

+  AxAyAzQ  (81) 
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(83) 


where  Catm  is  the  atmospheric  surface  pressure  (expressed  in  meters  of  water),  Ctp  is  the  tidal 
potential,  pt  is  the  internal  (baroclinic)  pressure,  ua  and  v°  are  the  horizontal  advection  velocities, 
and  F*  and  F*  are  the  finite  difference  forms  of  the  horizontal  momentum  mixing  terms.  The 
variables  are  evaluated  at  time  level  n  unless  otherwise  noted.  Note  that  the  temporal  changes  in 
the  height  of  the  grid  cells  Az  that  occur  on  the  sigma  coordinate  part  of  the  grid  are  accounted 
for  in  the  finite  difference  equations. 


The  primary  grid-cell  dimensions  Ax,  Ay,  and  Az  are  defined  at  the  center  of  the  grid  cells. 
The  superscripts  u,  v,  and  w  denote  grid-cell  dimensions  at  the  u,  v,  and  w  velocity  locations, 
respectively.  These  are  obtained  by  averaging  the  grid-cell  dimensions  of  the  adjoining  grid  cells 
(e.g.,  Axu  =  Axx). 

The  evaluation  of  the  surface  elevation  term  C*  can  be  distributed  among  any  of  the  three  time 
levels,  i.e., 


(*  =  aiCn+1+«2Cn+«3Cn_1, 


(84) 


where  the  temporal  weights  ai,  c*2,  and  are  specified  by  the  user  (see  Section  3.6). 
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The  variable  CCUrv  in  Eqs.  (79)  and  (80)  is  used  to  correct  the  horizontal  momentum  advection 
for  the  horizontal  curvature  of  the  grid.  The  curvature  correction  term  for  advection  has  a  form 
similar  to  that  of  the  Coriolis  term  (Appendix  A).  Ccurv  is  calculated  as 


_  -yMAy)  _  _x5j2y(Ax) 
2AxAy  2AxAy 


(85) 


On  the  sigma  coordinate  part  of  the  grid,  an  option  is  provided  to  calculate  the  horizontal 
diffusion  of  scalar  fields  relative  to  a  spatially  smooth  mean  or  climatological  field.  In  this  case,  T* 
and  S*  in  the  horizontal  diffusion  terms  are  set  equal  to  T  —  Tmean  and  S  —  Smean ,  respectively, 
where  Tmean  and  Smean  are  mean  or  climatological  or  horizontally  averaged  fields.  On  the  2-level 
grid  and  on  the  sigma  coordinate  grid,  if  this  option  is  not  used,  T*  =  T  and  S*  —  S. 

3.6  Calculation  of  the  Free-Surface  Mode 

The  free-surface  mode  is  calculated  implicitly.  Hence,  the  surface  pressure  gradients  and  the 
divergence  terms  in  the  surface  elevation  equation  have  a  component  at  the  new  time  level  being 
calculated.  The  finite  difference  equations  for  the  free-surface  mode  are 

A^^(^)  =  _AyuDug§x{aiCn+l  +  Q^n  +  ^n- 1)  +  (86) 

Ax^fyV  S2t(Dvv)  =  -  Axv  Dv  g5y(ai(n+1  +  a2(n  +  a3C"-1)  +  DVG~V,  (87) 

=  —Sx(  Ayu  (f3\(Duu)n+l  +  fo(Duu)n  + 

-  6y(Axv(f31(Dvv)n+1  +  fo(Dvv)n  +  fo(Dvv)n-1))  +  AxAyDQ,  (88) 

where  DUG Z  and  DVG „  are  the  vertical  integrals  of  all  the  terms  on  the  left  side  of  Eqs.  (79)  and 
(80),  respectively,  except  for  the  surface  elevation  gradient  terms,  and  Du  =  D  and  Dv  =  D  ' . 

The  variables  ai,  a2,  and  a3  are  constants  that  specify  the  fractional  weighting  of  the  surface 
elevation  gradient  in  the  momentum  equations  at  the  new,  current,  and  previous  time  levels.  Simi¬ 
larly,  /?i,  /?2,  and  /%  specify  the  fractional  weighting  of  the  divergence  terms  in  the  depth-averaged 
continuity  equation  at  the  new,  current,  and  previous  time  levels.  These  weightings  can  be  set  by 
the  user.  Commonly  used  values  are  ot\  =  a3  =  =  (33  =  0.5  and  0'2  =  do  =  0  (see  Section  4.4). 

The  equations  for  the  free-surface  mode  are  solved  by  substituting  the  expressions  for  ( Duu)n+l 
and  ( Dvv)n+1  from  Eqs.  (86)  and  (87)  into  Eq.  (88)  and  solving  for  the  new  surface  elevation  ("+1. 
The  resulting  equation  is  an  elliptic  equation  for  £n+1,  which  can  be  solved  with  an  iterative  or  a 
direct  method  (NCOM  currently  uses  an  iterative  solver). 

3.7  Baroclinic  Pressure  Gradient 

On  the  sigma  coordinate  part  of  the  grid,  the  baroclinic  pressure  gradient  is  calculated  for  a 
particular  layer  k  as 
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— SxPi\k  =  —SxPilk- 1  +  ~{]rDu{A(Tk- 1  +  +  Pk) 

Po  Po  Po  4 

-\(<rk-i  +  ak){6xD)(S^)\  (89) 

and  on  the  z-level  part  of  the  grid,  the  baroclinic  pressure  gradient  is  calculated  as 

— SxPi\k  =  —&xPi\k- 1  +  —\(kzk-i8xpk-\  +  A  zkSxpk),  (90) 

Po  po  Po  £ 

where  p*  =  p  —  p(z),  and  p(z)  is  the  horizontally  averaged  density. 

3.8  Horizontal  Advection 

The  advection  velocities  ua  and  va  in  Eqs.  (79)  to  (83)  are  calculated  from  un  and  vn,  respec¬ 
tively,  but  the  vertical  means  of  ua  and  va  are  adjusted  to  match  the  mean  vertical  velocity  needed 
to  account  for  the  change  in  the  surface  elevation  between  time  levels  n  —  1  and  n  +  1.  Hence, 

ua  =  un  +  f3lvn+l  +  02 un  +  -  un  (91) 


va  =  vn  +  (3ivn+1  +  02vn  +  (3zvn~l  -  vn.  (92) 

The  purpose  of  this  adjustment  is  to  ensure  that  the  advection  velocity  field  satisfies  continuity 
exactly  for  advection  of  the  scalar  fields. 


3.9  Horizontal  Mixing 


For  the  grid-cell  Reynolds  number  calculation  of  the  horizontal  mixing,  the  horizontal  mixing 
coefficients  are  calculated  at  the  staggered  velocity  points  as 


-  max 
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Avm  =  max 
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|t>n|Ayp' 
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and  the  finite  difference  form  of  the  mixing  terms  for  the  momentum  equations  is 
F*  =  5x{A^rA^r^JA^x  8xun~l)  +  8v(AxvAzvAlI/AyvX6yun-1), 


F*  =  8x{AyuAzuA\I/AxuU8xvn~l)  +  8y{AxvAzvAvM/Ay^8yvn-1). 


(94) 

(95) 

(96) 


For  the  Smagorinsky  horizontal  mixing,  the  horizontal  mixing  coefficients  are  calculated  at  the 
grid-cell  centers  as 

2  •  (97) 
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If  the  mixing  coefficients  are  then  averaged  to  the  staggered  velocity  points, 


A%[  —  Am  ? 

(98) 

ii 

(99) 

The  finite  difference  form  of  the  mixing  terms  for  the  momentum  equations  used  for  the 
Smagorinsky  scheme  is 

K  = 

+  5y{A^A^jA^X  5yun~l  +  AxvAzvAvMIAxvX8xvn~l),  (10°) 

F*  =  £x(Ay«AzM]{f/Ay«1'  5yun~l  +  A^A^A^/Aa^^”-1) 

+  ^(2A?MyA/<5/-1).  (101) 


3.10  Vertical  Mixing 

Vertical  mixing  is  fully  implicit,  e.g.,  the  vertical  heat  flux  is  computed  as 

firpn+l 

-KB-*—. 


(102) 


Fully  implicit  vertical  mixing  is  needed  to  avoid  spurious  flip-flopping  of  the  vertical  gradients, 
which  can  occur  with  a  partially  implicit  scheme  when  the  vertical  eddy  coefficients  become  very 
large  in  regions  of  strong  vertical  mixing.  The  implicit  treatment  of  vertical  mixing  couples  the  3-D 
prognostic  equations  in  the  vertical  and  requires  the  use  of  a  tridiagonal  solver  at  each  horizontal 
point  to  solve  for  the  new  values. 


The  finite  difference  form  for  the  prognostic  equations  for  the  MYL2.5  turbulence  scheme 
(Eqs.  (32)  and  (33))  is 
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The  variables  are  at  time  level  n  unless  otherwise  noted.  The  new  values  of  Kmx  and  Khx  for  the 
MYL2.5  scheme  are  then  calculated  as 


KMl  =(en+1qn+1SM  +  K%lId1)/2, 


(105) 


KHl  =  (£n+lqn+1SH  +  K%)/2,  (106) 

where  Kffi  and  are  the  values  of  KMl  and  Kh1  calculated  on  the  previous  timestep.  The 
purpose  of  averaging  with  values  calculated  on  the  previous  timestep  is  to  provide  a  strong  temporal 
smoothing  of  Kmi  and  Kux-  Without  this  type  of  averaging,  the  vertical  eddy  coefficients  tend  to 
be  noisy.  Su  and  Sh  are  calculated  from  G j;  using  Eqs.  (42)  and  (43)  where 


Gh  —  min 


0.028, 


/ pn+1 \ 

j2  9  dpn 

\qn+1, 

/  p0  dz 

(107) 


Note  that  the  eddy  coefficients  calculated  with  the  MYL2.5  scheme  in  NCOM  are  not  applied 
until  the  next  timestep  (Section  3.12).  Since  the  values  of  momentum  and  density  used  to  calculate 
the  new  values  of  KAIl  and  KHl  for  the  MYL2.5  scheme  are  at  time  level  n  and  the  new  values  of 
KMi  and  KHl  are  applied  at  the  next  timestep,  the  vertical  eddy  coefficients  are  calculated  from 
the  same  leapfrog  solution  (i.e.,  odd  or  even)  to  which  they  are  applied,  which  helps  to  suppress 
timesplitting. 

For  the  MYL2  scheme,  Sm ,  Sh,  and  q  are  calculated  using  Eqs.  (50)  to  (64)  with  the  values 
of  momentum  and  density  at  time  level  n  —  1.  For  the  MYL2  scheme  in  NCOM,  the  new  values 
of  Kmx  and  Kj!i  are  applied  on  the  same  timestep  at  which  they  are  calculated  (Section  3.12). 
Hence,  as  for  the  MYL2.5  scheme,  KhIl  and  KHl  are  calculated  from  the  same  leapfrog  solution 
to  which  they  are  applied.  Also  as  for  the  MYL2.5  scheme,  the  new  values  of  Kmi  and  Khx  are 
averaged  with  the  previously  calculated  values  to  reduce  noise. 

3.11  Bottom  Drag 

The  bottom  drag  is  partially  implicit  to  improve  stability  and  is  calculated  as 
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Kmjt-  =  cbun+1  |  v”-1  |  (108) 

oz 

and 

Km^  =  cbvn+1  |  v"-1  |  .  (109) 

OZ 

The  explicit  part  of  the  bottom  drag  terms  is  at  the  old  time  level  n  —  1  to  avoid  exciting  time¬ 
splitting. 

3.12  Calculation  Sequence 

The  calculation  sequence  for  the  model  is  as  follows: 

(1)  Advection  velocities  and  horizontal  eddy  coefficients  needed  for  momentum,  new  densities, 
and  new  baroclinic  pressure  gradients  are  calculated.  If  the  MYL2  turbulence  scheme  is  being  used, 
new  vertical  eddy  coefficients  are  calculated. 

(2)  New  3-D  horizontal  velocities  are  calculated  and  the  forcing  terms  from  the  3-D  momentum 
equations  are  vertically  integrated  to  provide  the  forcing  terms  needed  for  the  depth-averaged 
momentum  equations  that  are  used  to  calculate  the  free-surface  mode. 

(3)  The  depth-averaged  momentum  and  continuity  equations  are  solved  for  the  new  surface 
elevation  and  depth-averaged  velocities. 

(4)  The  new  3-D  velocity  fields  calculated  in  Step  2  are  corrected  by  adding  a  depth-independent 
adjustment,  so  that  their  vertical  mean  agrees  with  the  new  depth-averaged  velocities  calculated 
in  Step  3.  This  effectively  corrects  the  3-D  velocities  for  the  new  surface  elevation  gradient. 

(5)  The  velocity  field  that  will  be  used  to  advect  the  scalar  fields  is  calculated  by  adding  a 
depth-independent  adjustment  to  the  3-D  velocity  fields  at  time  level  n,  so  that  the  depth-averaged 
advection  velocities  are  consistent  with  the  depth-averaged  continuity  equation. 

(6)  New  values  of  the  scalar  fields  (T  and  S)  are  calculated  using  the  advection  fields  computed 
in  Step  4.  If  the  MYL2.5  turbulence  scheme  is  being  used,  the  turbulence  fields  are  updated  and 
new  vertical  eddy  coefficients  are  calculated.  These  fields  are  Asselin-filtered  as  the  new  values  are 
calculated. 

(7)  An  Asselin  filter  is  applied  to  the  velocity  and  surface  elevation  fields.  The  filtered  3-D 
velocities  are  then  corrected  to  be  consistent  with  the  filtered  depth-averaged  velocities  using  the 
same  procedure  as  in  Step  4. 

The  adjustment  of  the  advection  velocities  in  Step  5  ensures  that  the  velocity  field  used  to 
advect  the  scalar  fields  is  numerically  nondivergent.  This  is  necessary  to  avoid  spurious  sources 
and  sinks  when  using  the  flux  form  of  numerical  advection. 

It  is  desirable  that  the  velocity  field  used  for  advection  of  momentum  also  be  nondivergent 
and  consistent  with  the  change  in  surface  elevation.  However,  when  the  3-D  momentum  equations 
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and  the  forcing  for  the  free-surface  mode  are  calculated  in  Step  2,  the  new  surface  elevation  is  not 
yet  known.  Hence,  it  is  not  possible  for  the  momentum  advection  and  the  change  in  elevation  to 
be  fully  consistent  without  some  sort  of  iterative  process.  This  is  also  a  problem  for  free-surface 
models  such  as  POM  that  use  a  separate,  small  timestep  (i.e.,  the  split-explicit  scheme)  for  the 
free-surface  equations.  Iteration  of  the  solution  of  the  3-D  momentum  and  continuity  equations 
and  the  depth-averaged  equations  (Steps  1  to  5  above)  to  eliminate  this  inconsistency  is  provided 
for  in  NCOM.  However,  in  tests  that  have  been  conducted  to  date,  the  effect  of  iterating  to  remove 
the  slight  inconsistency  between  the  momentum  advection  and  the  change  in  surface  elevation  was 
not  significant. 

Another  difficulty  in  the  numerical  calculation  involves  the  partially  implicit  bottom  drag 
term.  If  the  bottom  drag  were  treated  explicitly,  the  implicit  vertical  mixing  and  the  bottom  drag 
calculation  would  be  numerically  decoupled  from  the  solution  of  the  depth-averaged  equations. 
However,  when  the  bottom  drag  calculation  involves  the  new  velocities,  the  bottom  drag  calculation 
and  the  solution  of  the  depth-averaged  equations  are  not  decoupled.  The  initial,  uncorrected 
estimate  of  the  new  3-D  velocities  will  be  involved  in  the  calculation  of  the  bottom  drag,  which  is 
part  of  the  forcing  term  for  the  free-surface  mode.  This  is  another  reason  why  it  is  important  that 
the  initial  calculation  of  the  new  3-D  velocities  (i.e.,  uncorrected  for  the  new  surface  elevation)  in 
Step  2  be  as  accurate  as  possible.  Currently,  in  the  initial  calculation  of  the  3-D  velocities,  the  new 
surface  elevation  is  estimated  from  the  horizontal  divergence  of  the  velocity  field  at  time  level  n. 

3.13  Shrinkwrapping  and  Slicing 

The  model  calculations  are  shrinkwrapped  in  the  x-direction  to  avoid  calculating  over  land 
areas  on  the  sides  of  the  domain.  What  this  means  is  that  the  calculations  start  at  the  first  sea 
point  and  end  at  the  last  sea  point  along  each  row  of  grid  cells  in  the  x-direction.  How  successful  this 
procedure  is  in  reducing  calculations  depends  on  the  particular  distribution  of  land  and  sea  areas 
in  the  model  domain.  For  many  coastal  problems,  shrinkwrapping  can  provide  a  useful  increase  in 
model  performance. 

Slicing  is  the  procedure  of  calculating  through  the  model  domain  in  x-z  slices  and  performing 
as  many  calculations  as  possible  on  each  x-z  slice  before  moving  on  to  the  next  one.  The  purpose 
of  this  is  to  reuse  variables  that  have  already  been  brought  into  high-speed  cache  memory  to  avoid 
the  slower  accesses  from  main  memory.  If  each  calculation  in  the  code  is  performed  over  the  entire 
domain,  variables  are  more  likely  to  be  flushed  out  of  high-speed  cache  before  they  are  reused. 
Each  timestep  of  the  model  requires  two  major  passes  through  the  model  domain,  one  to  update 
the  momentum  fields  and  another  to  update  the  scalar  fields.  How  much  the  x-z  slicing  improves 
program  speed  depends  on  the  relative  size  of  the  domain  and  the  particular  computer  being  used. 

4.  MODEL  LIMITATIONS 

When  using  a  model,  it  is  important  to  be  aware  of  its  limitations.  The  purpose  of  this  section 
is  to  discuss  some  of  the  limitations  of  the  physical  and  numerical  parameterizations  in  NCOM. 

4.1  Hydrostatic  Approximation 

Like  most  ocean  models,  NCOM  uses  the  hydrostatic  assumption  in  which  the  vertical  momen¬ 
tum  equation  is  reduced  to  a  balance  between  the  gravitational  acceleration  and  vertical  pressure 
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gradient  terms  (Eq.  (3)).  In  particular,  the  vertical  acceleration  terms  are  ignored.  This  is  a  fairly 
good  assumption  for  large-scale  oceanic  flows  where  the  vertical  acceleration  terms  are  usually  fairly 
small  relative  to  the  hydrostatic  terms.  The  assumption  is  less  valid  at  small  scales  where  vertical 
acceleration  terms  become  relatively  more  important. 

Casulli  and  Stelling  (1996)  illustrate  some  situations  where  the  hydrostatic  assumption  breaks 
down.  One  example  is  that  of  a  propagating  freshwater  plume.  As  the  plume  advances  into  the 
ambient  fluid,  there  is  a  convergence  near  the  surface  at  the  front  of  the  plume  where  the  water 
is  accelerated  downward  and  underneath  the  advancing  plume.  If  such  a  plume  is  modeled  with 
a  hydrostatic  model,  one  finds  that  as  the  horizontal  grid  resolution  is  increased,  the  downward 
velocity  at  the  front  of  the  plume  continues  to  increase  to  values  that  are  significantly  higher  than 
the  correct  values.  This  is  because  the  vertical  velocity  in  the  hydrostatic  model  is  calculated  strictly 
from  the  divergence  of  the  horizontal  flow  and  its  magnitude  is  limited  primarily  by  the  horizontal 
grid  resolution.  If  the  full  vertical  velocity  equation  were  being  used,  as  in  a  nonhydrostatic  model, 
the  other  terms  in  the  vertical  momentum  equation  would  act  to  limit  the  vertical  velocities. 

Another  example  is  that  of  propagating  internal  waves.  Large  amplitude  internal  waves  tend 
to  steepen  as  they  propagate  due  to  amplitude  dispersion.  Solitons  form  due  to  a  balance  between 
the  steepening  effect  of  amplitude  dispersion  and  the  countering  effect  of  frequency  dispersion  in 
which  shorter  waves  travel  more  slowly  than  longer  waves.  Hydrostatic  models  do  not  account  for 
frequency  dispersion;  hence,  they  cannot  simulate  the  formation  of  soliton  waves. 

Nonhydrostatic  phenomena  can  have  scales  as  large  as  several  km.  For  example,  internal  soliton 
waves  can  have  horizontal  wavelengths  of  a  couple  of  km.  But  this  is  not  to  say  that  hydrostatic 
models  cannot  be  used  to  simulate  flows  at  small  horizontal  scales,  just  that  certain  processes  in 
which  vertical  accelerations  are  important  may  not  be  correctly  simulated.  Such  processes  are 
frequently  secondary  in  importance  to  investigating  the  variability  of  the  horizontal  flow,  a  task  for 
which  the  hydrostatic  model  generally  provides  useful  results. 

4.2  Sigma  Vertical  Coordinates 

Sigma  coordinates  provide  several  advantages  in  modeling  coastal  regions  including  the  ability 
to  accurately  represent  the  bathymetry  (given  sufficient  horizontal  resolution),  a  smooth  treatment 
of  bottom-following  flows,  the  ability  to  provide  increased  resolution  and  a  fairly  consistent  grid 
in  the  bottom  boundary  layer,  and  increased  vertical  resolution  in  shallower  water  (which  may 
be  desirable  if  the  focus  of  the  modeling  is  nearshore  processes).  The  major  problems  with  using 
sigma  coordinates  stem  from  inadequate  resolution  of  steep  slopes.  In  such  situations,  the  horizontal 
pressure  gradient,  advection,  and  mixing  terms  can  all  generate  significant  truncation  errors.  Note 
that  “steep  slope”  here  is  a  relative  term  and  can  refer  to  a  ship  channel  in  a  bay  in  a  high-resolution 
simulation  as  well  as  a  continental  slope  or  a  seamount  in  a  larger  scale  simulation. 

As  noted  in  Section  3.2,  the  horizontal  pressure  gradient  between  two  points  in  sigma  coor¬ 
dinates  is  calculated  as  the  difference  between  two  terms,  one  being  the  pressure  gradient  along 
the  sigma  layer  and  the  other  being  the  vertical  pressure  gradient  due  to  any  difference  in  depth 
between  the  two  points.  For  a  steeply  sloping  sigma  layer,  the  vertical  pressure  difference  can  be 
quite  large  and  the  net  horizontal  pressure  gradient  becomes  the  difference  between  two  large  terms, 
a  situation  in  which  numerical  truncation  errors  can  become  significant.  Subtracting  off  a  mean 
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vertical  density  gradient,  as  is  usually  done  and  is  discussed  in  Section  3.2,  helps  considerably,  but 
truncation  errors  can  still  be  significant. 


Haney  (1991)  discussed  the  “hydrostatic  consistency”  problem  with  the  horizontal  pressure 
gradient  in  sigma  coordinates.  The  degree  of  this  problem  is  defined  by  the  number  of  sigma  layers 
that  are  crossed  when  traversing  between  adjacent  horizontal  points  at  a  particular  depth.  If  a 
number  of  layers  are  crossed,  the  pressure  gradient  that  is  calculated  may  not  properly  represent 
the  actual  density  changes  since  the  changes  are  not  being  resolved  in  the  calculation.  Haney  argued 
that  it  was  undesirable  to  cross  more  than  one  sigma  layer  boundary  between  adjacent  horizontal 
points  at  a  fixed  depth  (it  can  be  seen  that  this  restriction  tends  to  be  most  strongly  felt  in  the 
deeper  layers).  This  turns  out  to  be  a  very  strict  limitation  to  which  few  modelers  using  sigma 
coordinates  actually  adhere. 


A  common  measure  of  the  severity  of  the  bathymetry  changes  is  defined  as  the  change  in  depth 
between  two  adjacent  horizontal  points  relative  to  the  mean  depth,  i.e., 
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(110) 


Haney’s  hydrostatic  consistency  limitation  for  the  bottom  layer  of  a  sigma  coordinate  model  can 
be  expressed  approximately  as 
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where  Acq,  is  the  fractional  thickness  of  the  bottom  sigma  layer.  Hence,  =  r/Acq,.  According  to 
Haney  (1991),  we  need  Tb  less  than  about  1.  Sigma  coordinate  modelers  frequently  use  a  criterion 
of  r  <  1.5  or  more,  which  for  typical  values  of  A ab  of  0.2  or  less  is  a  significant  violation  of  Haney’s 
hydrostatic  consistency  criteria. 


Martin  et  al.  (1998)  discussed  problems  with  advection  and  diffusion  along  sloping  sigma  sur¬ 
faces  when  using  sigma  coordinates.  A  problem  that  can  occur  with  the  advection  term  in  regions 
of  steep  slopes  is  that  two  adjacent  points  within  a  sigma  layer  may  lie  on  different  sides  of  a  ther- 
mocline  or  halocline.  In  this  case,  the  temperature  or  salinity  field  appears  as  a  sharp  front  to  the 
horizontal  advection  term.  If  there  is  persistent  advection  across  this  front  using  the  second-order, 
centered  advection  scheme,  there  will  be  what  is  referred  to  as  an  “advective  overshoot”  on  the 
upstream  side  of  the  front.  Such  an  overshoot  can  reach  about  one  third  the  magnitude  of  the  jump 
across  the  front  and  can  result  in  the  model  calculation  becoming  unstable.  The  problem  is  that 
the  front  is  not  adequately  resolved.  The  second-order,  centered  advection  scheme  requires  a  front 
to  be  resolved  by  at  least  a  few  gridpoints  to  avoid  significant  advective  overshoots.  The  remedy 
in  this  case  is  to  increase  the  horizontal  grid  resolution  or  reduce  the  steepness  of  the  bottom  slope 
so  that  the  fronts  are  better  resolved. 


As  illustrated  by  Paul  (1994),  diffusion  along  sloping  sigma  layers  can  effectively  act  as  vertical 
diffusion.  Since  horizontal  mixing  in  ocean  models  is  typically  much  larger  than  vertical  mixing, 
the  strong  diffusion  along  a  sloping  sigma  layer  can  result  in  excessive  vertical  diffusion  and  can  ex¬ 
cessively  erode  vertical  gradients.  The  practice  of  calculating  the  diffusive  fluxes  using  the  anomaly 
based  on  some  mean  or  climatological  field  (Section  3.2)  significantly  reduces  this  spurious  vertical 
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diffusion.  However,  if  there  is  a  large,  local  anomaly  relative  to  the  mean  field  being  subtracted  off, 
significant  spurious  diffusion  can  still  occur.  For  example,  in  an  area  of  downwelling  near  a  coast 
there  will  be  a  warm  anomaly,  which  can  spuriously  warm  the  shallower  water  closer  to  the  coast 
by  diffusion  of  heat  along  the  sigma  layers.  This  problem  can  be  reduced  by  accounting  for  the 
anomaly  in  the  mean  field  being  used  as  a  reference  (which  may  require  updating  the  mean  field 
periodically),  by  increasing  the  horizontal  grid  resolution,  or  by  decreasing  the  bottom  slope. 


The  common  answer  in  dealing  with  sigma  coordinate  problems  over  steep  slopes  is  to  increase 
the  horizontal  grid  resolution  to  better  resolve  the  slope  or  to  (artificially)  reduce  the  severity  of 
the  slope.  As  noted  by  Haney  (1991)  and  others,  increasing  vertical  resolution  does  not  help  and 
can  even  make  the  problem  worse.  What  the  maximum  slopes  are  that  can  be  tolerated  in  sigma 
coordinates  is  not  clear,  probably  because  the  answer  is  somewhat  situation  dependent.  Modelers 
typically  modify  their  bottom  slopes  as  little  as  possible  when  trying  to  alleviate  apparent  sigma 
coordinate  problems,  so  that  they  maintain  the  most  accurate  possible  bathymetry.  However,  it 
is  not  always  easy  to  discern  whether  or  not  steep  slope  areas  are  causing  spurious  results.  In 
some  cases,  running  grid  convergence  tests  or  comparing  against  a  different  (e.g.,  2- level)  vertical 
coordinate  simulation  may  be  the  only  way  to  determine  if  there  is  a  problem. 


4.3  Z-level  Vertical  Coordinates 


The  z- level  grid  in  NCOM  truncates  the  bottom  depth  to  the  nearest  model  level.  This  is  the 
simplest  way  to  implement  a  2-level  grid  in  an  ocean  model  but,  with  such  a  treatment,  the  accuracy 
of  the  representation  of  the  bathymetry  in  the  model  depends  on  the  vertical  grid  resolution  that  is 
used.  Martin  et  al.  (1998)  discussed  some  of  the  problems  resulting  from  the  use  of  such  a  2-level 
grid.  Basically,  the  model  generates  the  solution  for  the  stepwise  bathymetry  that  is  used  in  the 
model  rather  than  for  the  actual  bathymetry.  This  may  sound  obvious,  but  modelers  tend  to  think 
in  terms  of  the  bathymetry  they  are  representing  rather  than  the  actual  bathymetry  being  used  in 
the  model. 


The  effects  of  a  stepwise  bathymetry  are  readily  noticeable  in  a  barotropic  onshore  or  offshore 
flow  (Martin  et  al.  1998).  The  convergence  of  the  horizontal  flow  occurs  at  the  faces  of  the  bathym¬ 
etry  steps,  which  produces  a  vertical  velocity  “jet”  at  the  faces  of  the  steps.  If  the  steps  were 
actually  there,  this  would  be  the  correct  solution.  However,  if  the  steps  represent  an  approximation 
to  a  smoothly  varying  bottom  depth,  the  vertical  velocity  field  will  not  vary  smoothly  in  the  hor¬ 
izontal  as  it  should.  Increased  vertical  resolution  helps  to  reduce  this  problem.  A  better  solution 
would  be  to  truncate  the  bottom  grid  cells  on  the  2-level  grid  to  the  bathymetry. 


Since  the  2-level  grid  follows  level  surfaces,  advection  and  diffusion  along  the  2- levels  are  di¬ 
rected  horizontally.  If  isopycnal  surfaces  depart  from  the  horizontal  plane,  horizontal  diffusion  along 
the  level  surfaces  can  result  in  excessive  cross-isopycnal  diffusion.  This  tends  to  be  less  of  a  problem 
than  diffusion  along  steeply  sloping  sigma  layers,  but  it  can  still  be  a  problem,  especially  on  longer 
timescales  as  the  effects  of  spurious  diffusion  increase.  Simulations  of  basin-scale  circulation  with 
2-level  models  have  had  trouble  with  this,  which  has  spurred  the  development  and  use  of  isopycnal 
coordinate  models  (Bleck  et  al.  1992). 
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4.4  Implicit  Treatment  of  the  Free  Surface 

The  implicit  treatment  of  surface  waves  used  in  NCOM  is  not  as  accurate  as  the  split-explicit 
scheme  used  by  POM.  This  is  mainly  because  of  the  much  larger  timestep  used  with  the  implicit 
scheme.  The  error  in  propagating  surface  waves  involves  both  damping  and  phase  speed  error 
(Martin  et  al.  1998). 

NCOM  provides  user-selectable  temporal  weightings  for  the  surface  elevation  gradients  in  the 
momentum  equations  and  the  divergence  terms  in  the  depth-averaged  continuity  equation,  which 
are  the  principal  terms  involved  in  surface  wave  propagation  (Section  3.6).  Dukowicz  and  Smith 
(1994)  discuss  selecting  weightings  for  these  terms  for  an  implicit  treatment  of  the  free  surface. 
Dukowicz  and  Smith  used  an  equal  weighting  at  the  three  time  levels  for  the  surface  elevation 
gradients  and  a  fully  implicit  weighting  (i.e.,  fully  at  time  level  n  +  1)  for  the  divergence  terms. 
This  was  based  on  a  compromise  between  accuracy  and  noise  damping  requirements  for  their  large- 
scale  simulations.  The  Estuarine  and  Coastal  Ocean  Model,  Semi  Implicit  (ECOM-si)  version 
(Blumberg  1992)  uses  a  fully  implicit  treatment  of  both  sets  of  terms;  however,  this  treatment 
results  in  strong  damping  of  the  surface  waves  (Martin  et  al.  1998).  The  NRL  Layered  Ocean 
Model  (NLOM)  (Wallcraft  1991)  uses  an  even  split  between  the  old  and  new  time  levels. 

The  minimum  damping  of  surface  waves  with  the  implicit  scheme  can  be  obtained  by  using 
the  same  temporal  weighting  at  the  old  and  new  time  levels  for  the  surface  pressure  gradient  and 
depth-averaged  divergence  terms  in  the  free-surface  Eqs.  (86)  to  (88).  With  this  weighting,  there  is 
no  damping  by  the  leapfrog  scheme  itself;  the  damping  is  due  only  to  the  Asselin  filter  (Martin  et  al. 
1998).  The  damping  from  the  Asselin  filter  can  be  reduced  by  reducing  the  Asselin  filter  coefficient; 
however,  the  filter  coefficient  must  be  kept  large  enough  to  suppress  timesplitting.  If  the  damping 
of  the  surface  waves  must  be  reduced  further  or  if  the  phase  speed  error  for  the  propagation  of  the 
surface  waves  must  be  decreased,  then  the  model  timestep  must  be  reduced. 

4.5  Second- Order  Centered  Advection 

As  already  mentioned,  the  second-order  centered  advection  scheme  suffers  from  advective  over¬ 
shoots  at  sharp  fronts.  This  scheme  assumes  that  the  field  being  advected  varies  smoothly  in  space 
and  that  significant  gradients  are  resolved  by  at  least  a  few  points.  The  advantages  of  this  scheme 
are  that  it  is  simple,  computationally  efficient,  and  gives  fairly  good  accuracy  when  the  gradients 
are  well-resolved.  This  scheme  has  certainly  been  used  more  than  any  other  in  ocean  models  (e.g., 
in  the  Bryan- Cox  model  (Bryan  1969),  POM  (Blumberg  and  Mellor  1987),  NLOM  (Wallcraft  1991), 
and  many  others). 

Other  problems  with  the  second-order  centered  advection  scheme  that  are  related  to  the  over¬ 
shooting  problem  are  that  it  does  not  maintain  the  monotonicity,  extrema,  or  positive  definiteness 
of  the  advected  field.  These  are  sometimes  desirable  properties,  e.g.,  salinity  and  biological  con¬ 
stituents  being  advected  by  an  ocean  model  should  not,  in  principle,  take  on  negative  values. 

The  second-order  centered  advection  scheme  usually  requires  a  certain  amount  of  horizontal 
diffusion  to  maintain  smooth  solutions  in  strongly  advective  flows.  This  need  is  addressed  fairly 
directly  by  the  grid-cell  Reynolds  number  horizontal  mixing  parameterization  (Section  2.4).  Using 
a  small  value  of  the  grid-cell  Reynolds  number  (e.g.,  10)  usually  maintains  a  smooth  solution  but 
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tends  to  significantly  damp  the  flow  and  dissipate  sharp  gradients.  Using  a  larger  value  (e.g., 
100)  results  in  much  less  damping  but  can  result  in  a  significant  amount  of  small-scale  noise.  The 
Smagorinsky  mixing  scheme  is  less  directly  pointed  at  filtering  numerical  noise,  but  as  a  practical 
matter,  the  extent  to  which  it  suppresses  noise  and  damps  the  flow  can  be  regulated  by  changing 
the  value  of  its  scaling  parameter  Csmag  (Section  2.4). 

There  are  different  philosophies  regarding  the  presence  of  small-scale  noise  in  a  numerical  model 
simulation.  One  philosophy  is  that  all  such  noise  is  objectionable.  Another  is  that  that  there  is 
a  fair  amount  of  noise  in  the  real  ocean  anyway  and,  as  long  as  the  dynamical  processes  being 
simulated  are  not  significantly  affected,  a  certain  amount  of  noise  is  tolerable. 

4.6  Timestep  Limitations 

A  basic  principle  for  the  stability  of  explicit  numerical  schemes  (at  least  those  being  used  in 
NCOM)  is  that  a  signal  (i.e.,  the  propagation  of  a  wave  or  the  advection  or  diffusion  of  a  field) 
cannot  travel  more  than  a  single  grid  interval  in  a  single  timestep.  The  processes  that  limit  the 
timestep  used  in  NCOM  are  internal  wave  propagation,  horizontal  and  vertical  advection,  and 
horizontal  mixing.  Since  the  surface  waves  and  vertical  mixing  are  treated  implicitly,  they  are  not 
subject  to  a  timestep  limitation  with  regard  to  numerical  stability. 

The  speed  of  internal  waves  depends  on  the  vertical  density  stratification  and  the  depth.  In  the 
deep  ocean,  the  maximum  speed  of  internal  waves  is  typically  2  to  3  m/s.  In  coastal  areas,  where 
the  water  is  relatively  shallow  and  the  stratification  tends  to  be  weaker,  the  maximum  internal 
wave  speeds  are  usually  less  than  1  m/s.  The  horizontal  advection  speeds  depend  on  the  situation 
(e.g.,  on  the  strength  of  the  tidal  currents  and  the  winds)  but  are  typically  less  than  2  m/s  and 
may  be  less  than  1  m/s  in  weakly  forced  situations.  In  general,  the  maximum  speeds  for  internal 
wave  propagation  and  advection  tend  to  be  comparable. 

Vertical  advection  is  another  matter.  If  one  is  trying  to  use  high  vertical  resolution,  vertical 
advection  can  be  the  limiting  factor  in  setting  the  timestep.  With  sigma  coordinates,  the  vertical 
grid  spacing  is  reduced  in  shallow  water  and  this  reduced  grid  spacing  may  limit  the  timestep.  The 
use  of  a  stretched  vertical  grid  with  higher  vertical  resolution  at  the  surface  and/or  bottom  (on 
a  sigma  coordinate  grid)  and  lower  resolution  at  mid  depth  helps  the  situation  since  the  vertical 
velocity  must  go  to  zero  at  the  free  surface  and  at  the  bottom.  Another  factor  that  helps  this 
situation  is  that  vertical  mixing  in  shallow  water  due  to  winds  and  bottom  drag  helps  suppress  the 
generation  of  vertical  motion  relative  to  the  sigma  surfaces. 

Horizontal  diffusion  is  almost  never  a  factor  in  setting  the  timestep  since  the  horizontal  diffusive 
time  scale  necessary  to  achieve  smooth  solutions  is  usually  much  smaller  than  the  advective  and 
internal  wave  timestep  restrictions. 

For  coastal  simulations,  a  conservative  estimate  of  the  maximum  internal  wave  and  advection 
speeds  is  generally  about  2  m/s.  With  2  m/s  as  a  maximum  velocity,  the  maximum  allowable 
timestep  for  NCOM  for  a  1-km  grid  would  be  about  250  s  (the  effective  timestep  for  the  leapfrog 
scheme  is  2 At).  For  a  larger  or  smaller  grid  spacing,  the  maximum  timestep  would  change  propor¬ 
tionally.  Note  that  these  numbers  just  provide  an  approximate  value.  Vertical  advective  processes 
may  require  a  smaller  timestep,  and  it  is  possible  for  wave  and  advective  speeds  to  combine  to 
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generate  larger  effective  signal  velocities  than  either  advection  or  wave  propagation  considered 
separately. 

Another  consideration  for  the  timestep,  besides  stability,  is  accuracy.  Because  the  advection, 
baroclinic  pressure  gradient,  and  Coriolis  terms  are  treated  explicitly  and  are  centered  in  time  with 
the  leapfrog  scheme,  temporal  truncation  errors  are  small  and  temporal  accuracy  is  quite  good  as 
long  as  the  timestep  is  sufficiently  small  that  the  numerical  scheme  remains  stable  (and  the  inertial 
period  is  well-resolved  for  the  Coriolis  term).  Since  surface  waves  and  vertical  mixing  are  treated 
implicitly,  their  temporal  accuracy  must  be  given  separate  consideration. 

The  accuracy  of  the  implicit  treatment  of  the  free  surface  is  discussed  in  Section  4.4.  The 
damping  of  surface  waves  can  be  minimized  for  a  given  timestep  by  appropriate  selection  of  model 
parameter  values.  Further  reduction  of  surface  wave  damping  and/ or  reduction  of  the  surface  wave 
phase  speed  error  requires  decreasing  the  timestep. 

The  fully  implicit  treatment  of  vertical  mixing  is  subject  to  temporal  truncation  error.  In 
most  conditions,  this  error  in  the  mixing  rate  is  small  and  is  not  especially  noticeable  since  vertical 
mixing  tends  to  be  fairly  rapid  relative  to  the  usual  timescales  of  concern.  (There  is  also  the  point 
that  the  rate  of  vertical  mixing  is  probably  not  simulated  very  accurately  anyway  due  to  the  lack 
of  a  physical  model  that  correctly  accounts  for  all  the  significant  processes  involved.)  However,  in 
cases  with  very  large  vertical  mixing  coefficients  (e.g.,  strong  winds)  and  fine  vertical  resolution, 
it  is  possible  to  have  sufficient  temporal  truncation  error  in  the  vertical  mixing  such  that  spurious 
mixing  occurs.  Such  spurious  mixing  effects  were  observed  in  some  idealized  forcing  experiments 
where  high  vertical  resolution  was  used  but  have  not  been  noticed  in  any  realistic  model  similations 
to  date.  It  is  planned  to  investigate  this  mixing  problem  more  thoroughly. 

4.7  Drying  of  Grid  Cells 

If  the  surface  elevation  reaches  the  bottom  of  a  grid  cell  or  the  bottom  of  the  sigma  coordinate 
grid  (za),  NCOM  will  “blow  up.”  This  is  sometimes  the  reason  why  an  otherwise  smoothly  running 
simulation  suddenly  stops.  To  keep  this  from  happening,  the  minimum  depth  of  the  model  grid  cells 
(or  za)  must  be  set  deeper  than  the  minimum  depth  of  the  free  surface  expected  at  that  location. 
Of  course,  one  doesn’t  always  know  what  the  minimum  depth  of  the  free  surface  will  be,  so  that 
an  educated  (conservative)  guess  may  be  needed  when  initially  setting  the  minimum  depths  for  a 
particular  simulation. 

5.  PLANS 

For  future  development  of  NCOM,  it  is  planned  to  address  some  of  the  limitations  discussed 
in  the  previous  section.  Some  enhancements  that  are  being  considered  are  as  follows: 

•  Provide  the  option  for  an  advection  scheme  that  provides  some  degree  of  upwinding  at  fronts 
to  avoid  or  reduce  the  overshooting  problems  of  the  second-order  centered  scheme.  A  number 
of  ocean  modelers  have  begun  using  advection  schemes  that  use  some  degree  of  upwinding 
at  fronts  to  avoid  or  reduce  advective  overshoots.  Many  schemes  have  been  discussed  in  the 
literature  (Rood  1987;  Pietrzak  1995),  and  the  different  schemes  have  particular  advantages  and 
disadvantages.  They  all  tend  to  involve  significantly  more  computation  than  the  second-order 
centered  scheme. 
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•  Truncate  the  bottom  grid  cell  to  the  bathymetry  on  the  0-level  grid  to  more  accurately  represent 
the  bottom  depth. 

•  Implement  a  flooding  and  drying  scheme.  This  would  provide  several  benefits  including  (a) 
preventing  the  accidental  drying  out  of  a  grid  cell  from  prematurely  terminating  a  model  run, 
(b)  providing  more  accurate  simulations  in  shallow  coastal  areas  by  allowing  depths  and  land-sea 
boundaries  to  be  set  without  being  restricted  by  the  need  to  avoid  grid  cells  drying  out,  and  (c) 
providing  the  ability  of  the  model  to  accurately  simulate  the  flooding/drying  process  itself  (if 
an  accurate  flooding/drying  scheme  based  on  correct  physics  is  developed). 

6.  SUMMARY 

This  report  provides  a  description  of  NCOM  Version  1.0.  The  model  has  a  free  surface  and  is 
based  on  the  primitive  equations  and  the  hydrostatic,  Boussinesq,  and  incompressible  approxima¬ 
tions.  A  choice  of  a  grid-cell  Reynolds  number  or  the  Smagorinsky  scheme  is  provided  for  horizontal 
mixing,  and  a  choice  of  the  MYL2  or  MYL2.5  turbulence  models  is  provided  for  the  parameteri¬ 
zation  of  vertical  mixing.  An  option  is  also  provided  to  include  the  vertical  mixing  enhancement 
scheme  of  Large  et  al.  (1994)  to  parameterize  unresolved  mixing  processes  occurring  at  near-critical 
Richardson  numbers.  The  inclusion  of  a  source  term  in  the  model  equations  simplifies  input  of  river 
and  runoff  inflows. 

The  model  uses  an  Arakawa  C  grid,  is  leapfrog  in  time  with  an  Asselin  filter  to  suppress 
timesplitting,  and  uses  second-order  centered  spatial  finite  differences.  The  propagation  of  surface 
waves  and  vertical  diffusion  are  treated  implicitly.  The  horizontal  grid  is  curvilinear.  The  vertical 
grid  uses  sigma  coordinates  for  the  upper  layers  and  0-level  (constant  depth)  coordinates  for  the 
lower  layers,  and  the  depth  at  which  the  model  changes  from  sigma  to  0-level  coordinates  can 
be  specified  by  the  user  by  setting  a  single  parameter.  The  combined  vertical  coordinate  system 
provides  some  flexibility  in  setting  up  the  vertical  grid  and  easily  allows  comparisons  to  be  made 
between  simulations  conducted  with  sigma  and  2- level  coordinates. 

The  model  is  based  on  fairly  well-tested  ocean  model  physics  and  numerics.  However,  there  are 
a  number  of  limitations  of  the  model.  Some  of  these  are  mentioned  here  (the  model’s  limitations 
are  discussed  in  more  detail  in  Section  4). 

Since  the  model  is  hydrostatic,  vertical  motions  on  small  horizontal  scales  may  not  be  properly 
described.  This  does  not  prevent  the  model  from  being  applied  with  high  horizontal  resolution  to 
look  at  the  structure  of  predominantly  horizontal  flows.  However,  nonhydrostatic  processes  that 
can  occur  in  these  situations  will  not  be  correctly  simulated. 

Sigma  coordinates  can  accurately  represent  the  changing  bottom  depth  but  can  suffer  from 
truncation  errors  in  their  horizontal  advection,  diffusion,  and  baroclinic  pressure  gradient  terms  if 
steep  bottom  slopes  are  not  adequately  resolved.  The  solution  to  this  problem  is  to  increase  the 
horizontal  grid  resolution  or  artificially  decrease  the  severity  of  the  slope. 

The  0-level  grid  does  not  suffer  from  these  problems  but  has  limitations  of  its  own.  Since  the 
0- level  grid  used  in  NCOM  truncates  the  bathymetry  to  the  nearest  0-level,  the  accuracy  of  the 
representation  of  the  bathymetry  on  the  0-level  grid  depends  on  the  vertical  grid  resolution,  and 
the  stepwise  structure  of  the  0-level  grid  can  cause  some  distortion  of  flows  that  cross  the  steps. 
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Also,  the  2-level  grid  does  not  provide  very  consistent  resolution  in  the  bottom  boundary  layer 
unless  a  large  number  of  levels  are  used  over  the  depth  range  at  which  the  bottom  boundary  layer 
exists. 

The  second-order  centered  advection  scheme  provides  fairly  good  accuracy  for  advection  of 
fields  in  which  the  gradients  are  moderately  well-resolved  but  can  generate  advective  overshoots 
at  sharp  fronts.  In  setting  the  timestep  for  the  model,  the  timestep  limitation  for  the  propagation 
of  internal  waves  and  for  horizontal  and  vertical  advection  must  not  be  exceeded  or  numerical 
instability  may  result.  Also,  the  drying  out  of  a  grid  cell  due  to  depression  of  the  free  surface  to 
the  sea  bottom  in  shallow  water  or  the  bottom  of  the  sigma  grid  can  cause  a  model  simulation  to 
suddenly  terminate. 
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Appendix  A 


EQUATIONS  IN  ORTHOGONAL,  CURVILINEAR, 
HORIZONTAL  COORDINATES 


The  derivation  of  vector  differential  expressions  in  orthogonal,  curvilinear  coordinates  is  dis¬ 
cussed  in  Batchelor  (1970).  The  form  of  the  basic  equations  in  orthogonal,  curvilinear,  horizontal 
coordinates  is  taken  from  Blumberg  and  Herring  (1987)  and  Song  and  Haidvogel  (1994).  The 
notation  used  here  is  similar  to  that  used  by  Batchelor  and  by  Blumberg  and  Herring. 

Let  the  new  horizontal  coordinates  be  f  ]  and  6  and  let  a  and  b  be  unit  vectors  paiallel  to 
£1  and  6>  respectively.  Then  the  change  in  the  horizontal  position  vector  x  corresponding  to 
increments  in  £1  and  £2  is 

4x  =  hi<5£ia  +  h25^2^>i  (Al) 

where  hi  and  h2  are  metric  coefficients  that  scale  the  distance  along  the  coordinate  axes  and  are 
functions  of  £1  and  6- 

The  horizontal  velocities  in  the  new  coordinate  system  are  defined  (for  convenience)  to  be  u 
and  v  and  are  directed  along  £1  and  £2,  respectively,  i.e., 

U  =  hl~dt ’ 

v  —  /i2~T~.  (A3) 

at 


The  basic,  3-D  Eqs.  (1)  to  (6)  in  the  new  horizontal  coordinate  system  become 
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These  equations  in  curvilinear  horizontal  coordinates  look  very  similar  to  the  original  Eqs.  (1)  to 
(6)  in  Cartesian  horizontal  coordinates.  The  only  additional  terms  that  appear  in  the  curvilinear 
equations  above  are  the  curvature  correction  terms  for  advection  in  the  momentum  equations,  i.e., 
the  term 
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h\h2 
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in  the  equation  for  u  and  the  term 


1  (  dh2 

Jhfo  V  56 
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(All) 


in  the  equation  for  v.  These  terms  account  for  the  interchange  of  u  and  v  momentum  for  advection 
along  a  curvilinear  grid.  Since  u  and  v  are  directed  along  the  horizontal  coordinates  6  and  6, 
momentum  advection  along  these  coordinates  will  effectively  result  in  conversion  of  momentum 
between  u  and  v  if  the  coordinates  have  non-zero  curvature.  Note  that  the  advective  curvature 
correction  terms  in  Eqs.  (A4)  and  (A5)  are  combined  with  the  Coriolis  terms  since  they  have  a 
similar  form. 


The  Smagorinsky  horizontal  mixing  terms  for  momentum,  given  by  Eqs.  (25)  and  (26),  when 
converted  to  curvilinear  coordinates,  have  the  form  (Blumberg  and  Herring  1987) 
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The  form  of  the  terms  for  the  grid-cell  Reynolds  number  horizontal  diffusion  in  curvilinear  coordi¬ 
nates  is  analogous.  The  terms  involving  §|f  and  §|-  in  Eqs.  (A4)  and  (A5)  are  curvature  correction 
terms  for  the  horizontal  mixing  of  momentum.  These  terms  are  assumed  to  be  small  and  are  not 
included  in  POM  (Mellor  1996)  or  NCOM. 


Note  that  when  converting  the  equations  to  finite  difference  form  with  second-order  centered 
differences,  we  set 
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If  x!  is  taken  to  be  the  arc  length  in  the  6  direction  and  /iiA6  =  Ax'  is  taken  to  be  the  grid 
spacing  in  the  6  direction  and  then  the  primes  on  x  are  dropped,  then 
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Similarly, 
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where  Ay  and  A yu  are  the  grid  spacing  in  the  6  direction  at  a  grid-cell  center  and  at  a  u-point, 
respectively.  In  deriving  finite  difference  expressions  such  as  in  Eq.  (A20),  use  is  made  of  the  fact 
that  A6  and  A£2  in  the  transformed  coordinate  system  are  constants  (i.e.,  variations  in  the  grid 
spacing  are  accounted  for  in  the  metric  coefficients  hi  and  h2). 
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Appendix  B 

EQUATIONS  IN  SIGMA  VERTICAL  COORDINATES 


The  basic,  3-D  model  equations  written  in  curvilinear  horizontal  coordinates  and  sigma  vertical 
coordinates  are  (Blumberg  and  Herring  1987) 
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where  and  F(,  are  the  horizontal  mixing  terms  for  u  and  v  momentum. 
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